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1 Introduction 


Theories of the known, which are described by different physical ideas, may be 
equivalent in all their predictions and hence scientifically indistinguishable. However, 
they are not psychologically identical when trying to move from that base into the 
unknown. For different views suggest different kinds of modifications which might be 
made and hence are not equivalent in the hypotheses one generates from them in 
one’s attempt to understand what is not yet understood. 


R. P. Feynman [1966] 


There is a write-up that provides the steps to turn an equation-of-state 
(EOS) into a mass-radius relation of exoplanets. It must be visualized [66]! 
1.1 Three Modes of Compression 
Compression can be classified by its timescale in the following way: 

e Supersonic: adiabatic irreversible (Rankine-Hugoniot) 

e Sonic: adiabatic reversible (isentropic) 

e Static: isothermal 

Shock wave, by definition, travels faster than the sonic wave in the matter. 
Thus, matter is compressed before the sonic signal of compression arrives. 

1.2 Three Modes of Energy Transport 
Energy transport can be classified by its physical mechanism into three modes: 


e Convection: large-scale flow due to internal heat source and gravitational 
buoyancy 


e Conduction: propagation of vibrational (thermal) energy of neighboring 
atoms 


e Radiation: absorption, transfer, reflection, and propagation of photons 


All of them occur on a planet, but in different parts, and in different ways. 
Convection often occurs throughout the bulk interior of a planet. Conduction 
often occurs within a thin boundary layer. Radiation often occurs on or near 
the planet surface by interaction with host stellar radiation. 


1.3 Hydrostatic Equilibrium 


A planetary body, to a good approximation, is a spherical body with concen- 
tric layering with continuous or discontinuous physical and chemical properties 
changing along the radius (or equivalently, depth) in the planet interior. Here, 
we first consider a planetary body with uniform chemical composition. Thus, we 
leave aside the chemical change for now, and focus on the physical changes oc- 
curring in a planet interior. The physical changes include (1) continuous change, 
such as the continuous increase of pressure, temperature, and density, along an 
isentrope, and (2) discontinuous or abrupt change, such as a phase transition 
from fluid to solid. We will consider both cases in our calculation as follows. 

The most important physical parameter which plays in the mass-radius rela- 
tion of any object, is its density (p). As often cited in high-school mathematics 
textbook, it was probably Archimedes who first discovered the mathematical 
relationship between the radius (r) of a sphere, and its surface area (A = 4772), 
and its volume (V = tp), Then, volume is turned into mass by multiplying 
with a physical parameter (p). A good example is liquid water under room tem- 
perature and sea-level pressure (1 bar) has an approximate density of 1 gram 
per cubic centimeter (g/cc). In contrast, the air under room temperature and 
sea-level pressure (1 bar) has an approximate density of 0.001 gram per cubic 
centimeter (g/cc), which is one-thousand times less. 

For self-gravitating object that is larger than a few hundred kilometers, p 
will not stay constant but will increase towards the center. Let’s focus on p in 
between 0.01 and 10 g/cc. Below 0.01 g/cc, Hydrogen can be well approximated 
as Ideal Gas (diatomic molecules H3): 


p(density in g/cc) 


P(pressure in bar) = -83.14-(T(Kelvin)) (1) 


p(molecular weight in g/mol) 
That means even for the coldest case of ~100 K, at ~ 100 bar, Hə can still 
be treated as ideal gas (p ~ 0.01 g/cc ). 
Above 10 g/cc, Hydrogen can be well approximated as a degenerate con- 
densed matter, with the following EOS using density p and Temperature T as 
independent variables: 


Pota = Po + Pin + Pei (2) 


Piotal is the total pressure corresponding to a given pair of (p,T), for a specific 
composition. Pim is the thermal correction from (lattice) vibration of atoms. 
Pa is the electronic contribution if (partially) ionized or (partially) metallic. 
In between 0.01 g/cc and 10 g/cc, there is no simple analytic expression of the 
EOS. This range is where different methods of Quantum Molecular Dynamics 


(QMD) Simulations play an important role. This range of density happens to 
coincide with the average bulk densities of many planets. For convenience of 
future calculations, we express pressure (P) as a function of two independent 
variables: density (p) and temperature (T). We will later explain the reasons 
behind this choice of independent variables. 


1.4 Differential Equations 


Now, for planetary interiors, we treat every variable as dependent on m (the 
mass enclosed by the sphere of radius r inside the planet). One can think of m 
as a natural coordinate. For example, if a planet cools and contracts over time, 
then r would decrease but mass m is conserved. Therefore, we write density 
at a particular point in the planet as p[m] instead of pfr]. Likewise, we write 
temperature at a particular point in the planet as T[m] instead of T[r]. Even 
radius r is considered a function of m, written as r|m]. This choice of dependent 
and independent variables will help us set up the differential equations. 

Now, let’s turn out attention back to EOS for a moment. Since the density 
typically ranges over three orders of magnitude or more inside a planet (from 
center to surface, so does temperature, and pressure), it is convenience to the 
express the EOS in logarithmic relation. That is, to take the logarithmic of 
every variable (T,p, and P), and then express the EOS as the relationship 
among lgT, lgp, and IgP, etc. Take logarithmic on these variables also helps 
non-dimensionalize them, and facilitate dimensional analysis. 

For example, the ideal gas law turns into a simple linear relation after this 
manipulation: 


lgP = lgp + lgT + const (3) 


Since for mass-radius relation, the crucial parameter is density p. We write 
the hydro-static equilibrium equation and the mass-conservation equation into 
the following forms, which are explicit in lgp but implicit in lgP and lgT. 

mass conservation: 

plm]: (r[m])? - r'[m] = 1.83843 (4) 
hydro-static equilibrium: 
D 
4 lgp[lgp[m],lgT 
. 110. Jerllgelm)igTim] — —115.06 - 
(rim) - 3 | m 6) 


The two constants on the righthand side are introduced for the simplicity of 
unit conversion so that the mass and radius are expressed in Earth units, and 
pressure is expressed in GPa (10° Pascal). 


(Mẹ = 5.9742 * 1074kg) 
4r - 108 - (Rẹ = 6.371 * 108m)’ 


1.83843 == (6) 


(G = 6.67428 « 10711) - (109) - (Ma = 5.9742 * 104kg)? 


115.06 == 
DN 4r - (Re = 6.371 x 106m)4 


(7) 


In order to solve this set of two equations, we need to supply two more rela- 
tions, namely, the (1) EOS expressed as the functional dependence of lgp|lgp, lgT], 
and (2) the relation between lgp and IgT. 

The former relation of lgp[lgp, lgT] is a generic material-specific EOS. 

The latter relation is a one-to-one correspondences between lgp and lgT. It 
is based on our assumption of the thermal structure in the planetary interior. 
For isothermal assumption, lgT is constant, which essentially turns lgp to be a 
single variable function only dependent on lgp. For isentropic assumption, the 
tangential slope of EE is a dimensionless quantity called Griineisen parameter 
y: a material property typically on the order of unity along an isentrope. In most 
cases, the isentropic assumption is preferred for planet deep interior, for there is 
often energy source in the planetary interior that drives an interior convection 
and maintains the temperature profile close to adiabatic (isentropic). However, 
there are some exceptions, such as the existence of thermal boundaries at the 
upper and/or lower boundary of each homogeneous convective layer. Within 
the thermal boundary, the heat transport mechanism switches to conduction 
which renders a steeper temperature gradient. 

In this manner, the boundary condition for integrating the differential equa- 
tions is set as (1) prescribing a starting central density pp at the center, where 
m = 0 and also r = 0; (2) prescribing a surface density p,, where the outward 
integration stops. 

This approach avoids iterative solving for any aimed quantities, but will 
produce a final total mass and radius once the integration is completed. It will 
simultaneously solve for the density, temperature, and pressure profiles along 
the radial direction of planet interior, thus providing a full picture of the 1-D 
planet model. 

Then, by continuously changing the central density pg, one generates a set of 
mass-radius points, from which, one can generate a smooth mass-radius curve 
corresponding to the specific EOS. 


1.5 Griineisen Parameter 


Griineisen parameter, denoted as lower-case y, is an important dimensionless 
parameter that relates to both the slopes of isentropes and melting phase- 
transitions in a T-p phase diagram. We reserve the capital-letter I" for a different 
but related dimensionless quantity upcoming in a later discussion. 

It is defined as the isochoric derivative of Pressure P respect to specific 
internal energy u, multiplied by the specific volume v = 1/p: 


OP 1 OP 
D Em 
v p 


differentiated along an isochor 


a) 
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This is especially convenient, such as in Mie-Grüneisen EOS, if Pressure P 
is expressed by independent variable v and u: 


P= P(v,u) (9) 


so that the full differential of P can be expressed as: 


OP OP 


On the other hand, as we more often encounter, the specific internal energy 
is expressed in its natural variables v, and s (specific entropy, or entropy per 
unit mass): 

u = ulv, s) (11) 


and according to the first law of thermodynamics, the full differential of u 
equals: 


du = —P -dv +T -ds (12) 


Through Maxwell relations, it can be shown that y relates to the slope of an 
isentrope in T-p space as: 


_ Əla T _ p fet 
Gee), tC), 
—— 


N 
differentiated along a reversible adiabat 


It characterizes the fractional temperature increase (T/T) along an isen- 
trope (constant s-path) due to fractional compression (dp/p). 
We define another dimensionless parameter capital-I as the : 


o OlnP 
r= å (14) 


differentiated along a reversible adiabat 


Their ratio y/T gives yet another dimensionless parameter: 


= OlnT One Oln P (15) 
= OlnP ~ (alp Olnp 
s s 
differentiated along a reversible adiabat 


IR 


To compare y with I and 7/T and emphasize their differences and relations: 


ey= (21) is the logarithmic-slope of T — p along an isentrope. We 


s 
call it the Griineisen parameter or the first adiabatic index. 


eT= (see ) is the logarithmic-slope of P — p along an isentrope. We 


call it the second adiabatic index. 


OlnP 


e ~T= (3) is the logarithmic-slope of T — P along an isentrope. We 
call it the third adiabatic index. 


In principle, both I and y and y/T vary along an isentrope, but rather 
slowly. Therefore, in certain approximation, they can be treated as constants, 
and this shall bring huge advantage in considering qualitative models, such as 
a polytropic model. 


1.6 Isentrope (Reversible Adiabat) 
Along an isentrope, since ds = 0, Eq. [12] simplifies to du = —P - dv. Then, 


— 
differentiated along a reversible adiabat 

where the subscript s means the derivative is taken under constant s. 

However, as often the case in literature, the functional dependence of specific 
energy u on specific entropy s is not explicitly given. Rather, specific entropy s, 
or specific internal energy u is given as a function of density p and temperature 
T. Then, isentropic path can be calculated solely by knowing the functional 
dependence of P = P(T, p) and u = u(T, p) via: 


e integrative scheme, see [69], Appendix A: method for calculating entropy. 


u(T, "e ike 
To= 2- S (ER) (17) 
0:P0 
P / T / 
z u(T, p) DE / dp! , = i P(To,p") _ dT" i ME Å so 
T u p To To T 


where f = f(T, p) is the specific Helmholtz free energy. 


e differential scheme, see [18], Equation 2 therein: 


E) 
(er) (dr) [fP Ep Å 
a dnp /, dlgp / , pu (a) p 
ð 
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where the subscript p means the derivative is taken under constant density. 
Note that the quantity (2) is dimensionless. p, when multiplies to u, 
converts it from energy per unit mass, into energy per unit volume, which 


is the same unit as force per unit area, that is, pressure. 


Here we concentrate on the differential scheme. The full differential incre- 
ment in P = P(T, p) or u = u(T, p) is separated by the effects of increment in 
Temperature T or increment in density p = + as: 


OP OP 
dP = (3) -dT + (3) -dp (19) 
P T 
Ou Ou 
m=(2) arr (2) a 
P T 


When we fix the specific volume v or density p, that is moving along an 
isochor, as in the definition of Grüneisen Parameter y, then, dv = 0, or equiva- 
lently, dp = 0, and we have: 


Then Grüneisen Parameter y can be re-written as: 


er) (rm) . 


(ower) á (2uor) 


It should be noted that pressure P increases faster along an isentrope com- 
pared to an isotherm for the same amount of increase in density p: 


OP) _ (aP) (oP) (ar a 
Op), pjr OF] NG, 

—"—" =- 

isothermal part due to y 


multiply density p across both sides of Eq. we have: 


(srs), > (aims), er), 


— —— —S 


Ks, isentropic bulk modulus Kr, isothermal bulk modulus due to y 


Then, plugging Eq. [13] and Eq. [21] for Griineisen Parameter y, we have: 


always positive 


Therefore, K, > Kr and it is more difficult to compress matter along an 
isentrope, compared to an isotherm. 


1.7 Melting 


Melting, or fusion, refers to the phase transition from solid to liquid. The same 
process in the opposite direction from liquid to solid is called Solidification: 
Solid ==" — Liquid (25) 
Solidification 

On the microscopic scale, melting is associated with the breaking of cer- 
tain bonds in a crystal (Fig. o, and is sometimes described by the dislocation 
theory [77]. 

For a single-component system, melting is often described by a melting curve, 
in the P-T phase diagram or p-T phase diagram. 

When melting occurs, not only the Temperature T is fixed at Tm, but also 
the Pressure P is fixed at Pm. Therefore, the melting process is simultaneously a 
reversible isothermal and reversible isobaric process. However, there is typically 
a volume (density) change associated, and heat absorption (or release) involved. 

The differential relation between Py and Tm if we move along the melting 
curve is expressed by the Clapyron-Clausius relation as: 


ð Pm As Ah 
= = 2 
& i Av Tu -Av ( 6) 


along melting 


At high pressure, we are mostly concerned with the volume-driven fluid-to- 
solid phase transition due to compression, not the entropy-driven fluid-to-solid 
phase transition due to cooling. 

With sufficient amount of compression, a fluid transitions into a solid accom- 
panied by a finite (specific) volume change (Av). Unlike the gas-liquid phase 
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Figure 1: Different types of bonding in crystals. When melting occurs, some 
bondings are disrupted or dislocated, to allow flow. Copyright 2010 by Ency- 
clopædia Britannica, Inc. [2]. 


boundary that vanishes above the Critical Point C.P., fluid-solid phase bound- 
ary exists up to ultra-high pressure and temperature. The (specific) volume 
change (Av) means discrete density p jump along the melting curve. Therefore, 
viewed in the T-p diagram, the melting curve is actually represented by two 
separate curves, one representing the density of solid pgojiq along the melting 
curve, and the other representing the density of liquid piquia along the melting 
curve. The area in between these two curves are forbidden. 

argues that the slope of melting curve in the p-T phase diagram at 
high-pressure is related to Griineisen Parameter y as: 


OlnT 
(Sine)? em 
— 


— 
along melting 

where the subscript M means the path of differentiation is taken along the 
melting curve, that is, maintaining the phase equilibrium. 

In order to understand the melting process in more physical detail, let’s 
introduce the (specific) enthalpy Ah of melting (sometimes also known as en- 
thalpy of fusion), that is, the total energy required to melt a unit mass of matter 
at a given pressure Pm: 


Ah= hiquid KE hsolid (28) 


From the perspective of specific internal energy u and specific volume v, it 
can be written as: 


Ah = (u+ Puo) — (u+ Puo) (29) 
liquid solid 
= (inaua + tone) + Py - Ga = usoa ) 
= Au + Py : Av 


Au generally contains two parts: the increase in the kinetic energy of atomic 
motion Aux, and the increase in the potential energy of the atomic interactions 
among neighbors Aup. 

so that, 


Au = Aux + Aup (30) 

[78] attributes melting to the free proliferation of dislocations. Because melt- 
ing involves no temperature change T = Tm, then, Aux = 0, and, 

Au = Aup (31) 


From the perspective of specific entropy s, it can be written as: 


Ah = (ru å s) + (£u Å s) (32) 
liquid solid 


= Tum: (sii = st) 


= Tu - As 


argues that at high pressure, both solid and liquid have close-packing 
atomic arrangements. The difference between them is well-characterized by the 
dislocation theory. So the (specific) entropy increase of melting is: 


As = Aup/Tm (33) 


The interpretation of this formula is that the (specific) entropy increase As 
at melting is primarily contributed by the increase of spatial degree of freedom. 
The thermal degree of freedom stays roughly the same. Since for vibrations 
above a certain frequency (Frenkel frequency), molecules in a liquid behave 
like those in a solid, and can support shear phonons [26]. Thus, the heat ca- 
pacities are similar between the solid and liquid phases along the melting curve, 
provided that the higher frequencies of the thermal vibrations dominate. 
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However, for the melting of superionic H20 ice, because its protons are 
already mobilized and can flow freely like in liquid, As is less than what is 
expected for melting of a complete solid phase. Melting there only regards the 
oxygen atoms in a crystal lattice. Therefore, there is an additional factor of 3 
when calculating the slope of melting curve of the superionic H>0 ice in the 
p-T phase diagram at high-pressure: 


(mt) Gee ay 


Olnp 
—— 
along melting of superionic ice 


Please note that Eq. and Eq. are only approximate, and should be 
checked against experimental data or DFT-MD (Density Functional Theory- 
Molecular Dynamics) calculations. 


1.8 Ionization 


It is commonly known that under ambient conditions (25°C, 1 bar), the H2O 
molecules in the liquid dissociates into a small amounts H* and OH ions in a 
dynamical equilibrium: 


H20 Dissociation H+ i OH- (35) 
Recombination 


The degree of ionization, characterized by the mole fractions of ions, in- 
creases with both increased pressure and increased temperature : 


Mole fraction 


P (GPa) 


Figure 2: The concentrations of chemical species as predicted by the first princi- 
ples simulations in [44] are shown as a function of pressure. Inset: The lifetime 
of the H2O molecule is shown as a function of pressure. All other species had 
lifetimes between 30 and 50 fs, independent of pressure. Copyright 2009 by AIP 
Publishing. 
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Experiments show an enormous increase in Ah (enthalpy of fusion) along the 
melting line of H20 at around ~40 GPa [44]. They attribute this measurement 
to significant molecular dissociation of H2O upon melting (see Fig. |2). Thus, 
the fluid in this regime is termed tonic fluid. The tonic fluid will have higher 
electric conductivity (o in S/cm) compared to molecular fluid, owing to the 
mobility of ions. 

We expect similar phenomenon of increased molecular dissociation into ions 
to occur in fluid Ammonia (NH3) under increased pressure and temperature: 


NH; Dissociation Ht + NH2- (36) 


Recombination 


A schematic phase diagram of water and ammonia is shown as follows: 
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Figure 3: Schematic phase diagram for water and ammonia, or their mixture. 
Also shown are trajectories of isotherm, isentropes, and single/multiple-shock 
Hugoniot. A single-shock Hugoniot will likely miss the super-ionic regime. The 
detailed differences among the various ice structures/phases have been ignored 
as they are all approaching close-packing at sufficient pressures. 


12 


1.9 Electric Conductivity 
The electric conductivity å (in S/cm) is defined as: 


current density 


SN 


J 
E 
Ww 
electric field magnitude 
It can be caused by: 
e de-localized (mobile) ions 
e de-localized (mobile) protons (Ht) 


e de-localized (mobile) electrons (e7) 


The range of o spans over many orders-of-magnitude: 


resistivity e(Q-cm) 
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10-1810-1610-14 10-1210-1010-8 10-6 10-4 10-2 1 102 104 106 108 
conductivity o (S/cm) 
— insulator — —— semiconductor ——|+conductor — 
© 2004 Encyclopædia Britannica, Inc. 


Figure 4: Range of electric conductivity (ø in S/cm). Copyright 2004 by Ency- 
clopædia Britannica, Inc. [3] 


Here we concentrate on the electric conductivity resulting from electrons 
(67). 
Quote from the Encyclopædia Britannica [4] [B] [I]: 


Materials are classified as conductors, insulators, or semiconductors 
according to their electric conductivity. The classifications can be 
understood in atomic terms. Electrons in an atom can have only 
certain well-defined energies, and, depending on their energies, the 
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electrons are said to occupy particular energy levels. In a typical 
atom with many electrons, the lower energy levels are filled, each 
with the number of electrons allowed by a quantum mechanical rule 
known as the Pauli exclusion principle. Depending on the element, 
the highest energy level to have electrons may or may not be com- 
pletely full. If two atoms of some element are brought close enough 
together so that they interact, the two-atom system has two closely 
spaced levels for each level of the single atom. If 10 atoms interact, 
the 10-atom system will have a cluster of 10 levels corresponding to 
each single level of an individual atom. In a condensed matter, the 
number of atoms and hence the number of levels is extremely large; 
most of the higher energy levels overlap in a continuous fashion ex- 
cept for certain energies in which there are no levels at all. Energy 
regions with levels are called energy bands, and regions that have no 
levels are referred to as band gaps. 


From the perspective of energy, here we show a schematic drawing of the 
energy levels (band and gap), for different types of conductors, insulators, or 
semiconductors: Fig [5] 


Band Theon. C ig noring D-E intevartan thus, in ) 
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G > 
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metal — Semirmatal com tdu ole * sekk 


Figure 5: Band Theory. The envelope of bands can be understood as a repre- 
sentation of the density of states. The topology of band determines the nature 
of the type of electrical conductivity in the material. Here we only consider 
Pauli exclusion principle and have ignored e~-e~ interactions/coupling, thus, 
posing an incomplete theoretical picture. The probability of the occupancy of 
the energy levels fill up roughly to the Fermi energy Er, and can spill over a 
few kgT upwards. 
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In general, the temperature-dependence of o is different for a metal and an 
insulator (semi-conductor): 
For metal: 


T 
F= 69+ 7 XT (38) 
0 


For insulator (semi-conductor): 


Ea 
TE er) (39) 


where EF, is the activation energy for a given insulator (semi-conductor). The 
o of metals depends on the temperature and the ambient conditions much less 
because in metals, the state-of-conduction is a non-excited state. Meanwhile, 
in order to make an insulator (semi-conductor) conductive, it is necessary to 
transform them into an excited (activated) state [56]. 

In the next sections, we will show that both increase in Temperature and 
Density (or Pressure) will increase the conductivity of a material. See also 


in [84]. 


1.10 Metallization 
due to Density Increase 


1.10.1 Scaling Relations 


At high enough density, that is, with significant compression, everything is 
expected to become metal. The fundamental reason is simple: as atom spacing 
becomes small, their outermost electronic orbitals are forced to overlap, and 
electrons become de-localized and are shared with other atoms [47]. 

In terms of the band theory, the topology of band (Fig B) changes when the 
spacings between neighboring atoms change. 

At T = 0 Kelvin, the transition from non-metal to metal state can be esti- 
mated by the Mott’s Formula [39] [37]: 


n!/3- at œ 0.25 (40) 


Here ne the critical atom density that corresponds to a critical distance (de) 
between close-packing atoms, at which a first-order, discontinuous transition 
from a metal to a non-metal would occur at T = 0 K. ne has the dimension 
of reciprocal-volume, thus, its product with a distance a7 is dimensionless. af 
is an equivalent Bohr (hydrogenic) radius of the isolated (low density) atomic 
state (in this instance). a7; can be estimated as follows [39]: 


Ka Rh 

ag å = (41) 
mše 

where Ks is the static (low-frequency) dielectric constant of the host mate- 


rial, and mf! is the effective mass of an electron in the host conduction band. 
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At T > 0 Kelvin, say 1000 K, the scaling parameter ni! : 


from non-metal to metal state is relaxed to [37]: 


- az, for transition 


nt/3 . at, = 0.38 (42) 


Eq. and Eq. are convenient when we apply them to understand the 
non-metal-to-metal transition in the T-p phase diagram. For example, an eight- 
fold compression would double d. and increase both p and ne by eight-fold. 

Experiments and calculations show that almost all elements (not their com- 
pounds) become metal under the pressures of 200~300 GPa or above with the 
only exceptions of Helium, Carbon, Fluorine, Chlorine, and other noble gases 
or alike [47]. For example, metallic Oxygen is obtained in lab around 100 
GPa [38]. 


1.10.2 Fermi-Dirac Distribution for Electrons 


The electrons in a condensed phase can be divided into two groups; those that 
form the closed electron shells of the constituent atoms (the core electrons), 
and the remaining electrons of higher energy (the valence electrons). In metals 
the valence electrons can move more or less freely through the lattice with an 
effective mass which takes into account the influence of the periodic potential 
resulting from the lattice ions. One usually refers to them as the conduction 
electrons [45]. 

In the Sommerfeld model of a metal [76] in which the conduction electrons 
are assumed to form a gas of fermions with energies E and a density of states 
N(E). 

The probability of a state with the energy E being occupied by an electron 
in a state of thermodynamic equilibrium follows the Fermi-Dirac distribution 
function as: 


fo(E,T) = (ov EE = (or C= 1) (43) 


Here, Er is the Fermi level, defined as the Gibbs thermodynamic potential 
per particle, i.e., the chemical potential p at T = 0 Kelvin: Ep = p(0). Er 
numerically equals to the work that must be spent to add one particle (here 
one electron) to the system [56] [57]. Here we have ignored electron-electron 
interactions such as e~-e~ pairing which leads to super-conductivity. 

fo(E, T) experiences essential changes when the energy E is changed from 
~ (Ep — 2kpT) to ~ (Er + 2kpT). fo(E,T) can be integrated over the full 
momentum-space to find the mean energy (per electron) of all the electrons of 


a degenerate system: 
oo 2 2 
P p 2 
—: —,T | + (4 d 
f 2Me fo( ) Gu 
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In the two limiting cases, 
3kpT hen kpT E 
(E) 3 2 B when Kgl > BF, (45) 
sEr when kgT < Er. 


As can be seen from Eq. the energy of degenerate e gas at low tem- 
perature is independent of temperature, and thus, it does not contribute to 
the specific heat of the body, and this fully explains the Dulong-Petit law for 
solids [56]. (E) = 2 Er typically amounts to several electron-volts (eV). How- 
ever, at elevated temperature, when kgT > Er, the energy of the electron gas 
is determined by its temperature. The cross-over from quantum to classical 


statistics is set as follows: 


d= (rT) x1 (46) 


In other words, the dimensionless degeneracy parameter 0 = T/Tr is on the 
order of unity. Tp is the Fermi temperature defined from Fermi energy Ep as 
follows: 


E R? a 
= d I . Le. 
Tp = e ane de | T ") (47) 


where n is the number density of (conduction) electrons, that is, the number 
density of electrons in the conduction band, and m is the effective mass of an 
electron in the host’s conduction band. The Fermi energy of a metal is typically 
of the order of 2~10 eV [16] [15]. 

The matter in this cross-over regime is sometimes referred to as the Warm 
Dense Matter (WDM). It is the current interest of high-pressure shock wave 
experiments [68]. This cross-over regime is characterized by atomic matter 
for which potential energy of electron-ion interactions is comparable to kinetic 
energy of electrons. 


1.11 Plasma Phase Transition (PPT) 
due to Temperature Increase 


At high enough temperature, everything will become plasma. Plasma is a 
gaseous or gas-like phase characterized by free electrons (e7) due to thermal 
excitations. Owing to the free electrons (e7), it has high electric conductivity 
(o in S/cm) compared to ionic fluid. 

Plasma Phase Transition (PPT) is related with the metallization of fluid, as 
the PPT boundary sometimes continues smoothly into the denser regime in the 
Temperature-Density (T-p) phase diagram which then belongs to the liquid-like 
phase rather than a typical gas-like phase. The following is a simple schematic 
phase diagram for Hydrogen in T-p space: 

As seen in Fig. [6] the PPT draws the boundary between fluid molecular 
hydrogen (H2) and fluid metallic hydrogen (H* and free electrons). 
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Figure 6: Schematic Phase Diagram for Hydrogen. 


The fluid metallic hydrogen is immiscible with fluid helium (which is a noble 
gas) at low temperature, but eventually becomes miscible with fluid helium at 
higher temperature. 

According to calculations in upcoming chapter, Fig. (9 Fig. Fig. |13} 
and Fig. PPT of hydrogen occurs around 100 GPa, or correspondingly, 
at a density slightly below 1 g/cc. Jupiter’s interior may still be too hot to 
miss this miscibility gap. Pressures in Uranus’ and Neptune’s envelopes are 
everywhere lower than PPT at around 10° bar or 100 GPa. Therefore, PPT 
is most relevant to Saturn, and may explain its abnormal intrinsic luminosity 
and cooling history [79]. In more detail, PPT in fluid hydrogen also has a 
temperature-dependence, in addition to its dependence on Pressure P, as shown 
in the P-T phase diagram of dense hydrogen (Fig. |7| experiments show a rapid 
increase in electrical conductivity (o in S/cm) across PPT [67]). 


1.12 Mixing and Un-Mixing 
1.12.1 Thermodynamics of Mixing 


The fundamental reason for mixing and un-mixing is the differences in inter- 
molecular forces or specific molecular effects between different species, even 
though they are chemically non-reacting. 

Quantitatively, one can estimate the Gibbs free energy of mixing (AG mixing): 
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compression [68]. Copyright 2016 by American Physics Society. 


AG mixing => Gafter mixing T Gbefore mixing 


For a given Temperature T, and the number fraction of each individual 


= A Fixing =T. AS mixing 


= (Has mixing T Hbefore T =T. ASmixing 


component participating in the mixing, 


e If AGmixing < 0, then mixing is favored 


e If AGmixing > 0, then un-mixing is favored 


AH mixing is the enthalpy of mixing, and ASmixing is the entropy of mixing. 


AFmixing can be calculated from ab initio simulations : 
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A Fixing = Faster mixing (21 T25 085 ox) = Hbetore mixing (49) 


N 
= Hafter mixing CES -X ri H; 
t=1 
=AU +P. AV 


volume effect comes in here 


where the enthalpy of the mixture (Hafter mixing) is a function of £1, £2, ..., EN, 
which are the number fractions of every individual components participating in 
the mixing, assuming no chemical reactions occurring among them, and satis- 
fying the normalization condition: > = 1, 

Without a priori knowledge of the spatial distribution or orientation within 
the mixture, such as clustering or non-uniform distribution, ASmixing can often 


be approximated by the entropy of ideal mixing: 


ASmixing x ASideal mixing (50) 


N 
~—kp- > zi- ln zi 
i=1 


ASmixing is always positive, and it is weighted by Temperature T in Eq. 
Therefore, at high-enough temperature, mixing is always favored. AHmixing 
primarily depends on the energy state of the system. If the components interact 
with each other attractively, that is, by forming bonds and sharing electrons such 
as in metal, and then the total energy of the system is reduced after mixing, and 
AH mixing is small or even negative. On the other hand, if the components inter- 
act repulsively, then bringing them together cost a lot of energy and AH nixing 
is positive and large. 

One can think of Temperature T as playing the role of fulcrum on a lever, 
and the load on one side is AH, and the load on the other side is AS. We call 
this lever AG: 


AG=AH- T -AS (51) 
sowa 
fulcrum 
1.12.2 Relation to Metallization 


In particular, non-metal-to-metal transition plays a key role for mixing and un- 
mixing of different components in planetary interior, i.e., differentiation. The 
binary mixture of hydrogen-helium under pressure is just one example. Gener- 
ally speaking, metal (fluid) and metal (fluid) can readily mix with each other, 
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while metal (fluid) and non-metal (fluid) cannot. So mixing/un-mixing is viewed 
as the same problem as metallization. 

The story of compounds may be more complicated. For example, in silicates, 
the oxygen anions (having acquired e~s and silicon cations (having donated 
e7s), having assumed closed-shell structure and somewhat equivalent to noble 
gases, are much more stable than their elemental forms, and thus can remain 
non-metallic up to higher pressure (a few TPa or 1000’s GPa level). This is 
confirmed by both static-compression and shock-wave experiments. 

According to our earlier work on the simple scaling laws of planet interior 
pressure with its mass and material property [86], we find out that: 


e For rocky planets and water worlds, their characteristic interior pressure 
scales approximately linear in planet mass: 


2 
P~(M/Mo): (1- A (Ice mass fraction)) - 100 GPa (52) 


e For degenerate gas giants, their characteristic interior pressure scales ap- 
proximately quadratic in planet mass: 


P ~ (M/Myupiter)? © 1 TPa (53) 


For planets more massive than Earth, their core-mantle boundary of any 
sort may gradually become less distinct and more fuzzier with increasing mass, 
due to the metallization and increased miscibility of materials in their deep 
interior. The trends of mutual mixing of Fe-metal, silicates, ices, and all other 
materials dissolved in them increase. A planet core may still exist, because 
denser material is favored in terms of long-range gravitational energy towards 
the center. However, on the short-range, in terms of chemical gradation and 
entropy, mixing is favored more and more. 

This is the most natural explanation for Jupiter’s fuzzy core [83], as advo- 
cated by B. Militzer and D. Stevenson. In conclusion, the planets more massive 
than Earth will have different interior differentiation. The exact result is also 
evolution-history-dependent or path-dependent. 


1.13 Transport Properties 


Transport properties describe how a particular physical quantity is transferred 
through a system, for example diffusivity, viscosity, thermal conductivity, and 
electric conductivity (which we have just discussed). In recent years, there is 
a growing interest in calculating and measuring transport properties in high- 
pressure experiments, because they play important role in modelling and un- 
derstanding the physics of planetary interiors. 


1.13.1 Prandtl Number 


One important dimensionless parameter that characterizes the transport prop- 
erties is the Prandtl number Pr: 
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orderly motions dynamic viscosity density 


r — AN 
momentum diffusivity m lp 
Pr= = 
thermal diffusivity k /( c “p) 
NY Sø 
dis-orderly motions thermal conductivity 


specific heat capacity 
(54) 

It characterizes the relative importance between two energy transport modes: 
conduction versus convection. It then determines the relative thickness of a ther- 
mal boundary layer (dthermal boundary layer) Versus thickness of velocity boundary 
layer (dvelocity boundary layer) in a convective cell. 

Since there is no length-scale involved in Eq. 54 Pr is an EOS property which 
only depends on the type of fluid and the state the fluid is in. For most gases 
over a wide Temperature-Pressure range, Pr is approximately constant, due to 
the micro-physics involved that in the gas-phase both processes are caused by 
the kinetic motions of gas molecules. 


e Pr < 1 (Liquid Metals, Earth’s Core): Heat diffuses very efficiently and 
quickly. Thus, heat conduction is more significant compared to convection. 
Then, 


Öthermal boundary layer < Övelocity boundary layer (55) 


e Pr < 1 (Oils, Earth's Mantle): Heat diffuses very inefficiently and slowly. 
Thus, convection dominated. Then, 


Öthermal boundary layer > Övelocity boundary layer (56) 


e Pr ~ 1 (Many gases): Heat diffusion and convection are of comparative 
importance. Therefore, both effects need to be dealt with simultaneously, 
such as in aerodynamics. 


1.13.2 Rossby Number 


Another important dimensionless parameter that characterizes the large-scale 
transport properties is the Rossby number Ro. It characterizes the large-scale 
flow in the surface-layer of a planet, such as atmospheric motions and flow in 
ocean. 


characteristic length scale 


[Inertial Force lv Vo] U2/ aoa 
Ro = = ny (57) 
[Coriolis Force] |Q x v] U- Q 


angular frequency of rotation 


On Earth, Ro accounts for phenomena such as jet streams, long waves in 
westerlies, and polar fronts, etc. On Hot Jupiters, it shall account for the flow 
structure in and out of the hot-spot caused by host-stellar irradiation. 
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2 Variational Principle and Virial Theorem 


2.1 Basics: Hydrostatic Equilibrium Re-visited 


Owing to the simplicity of the form of specific internal energy differential along 
an isentrope (reversible adiabat) with constant specific entropy s = sọ: 


du = —P dv (58) 


We can construct a Lagrangian, which by means of variational principle, 
is equivalent to the differential equation of hydro-static equilibrium and mass- 
conservation [85]. 

Assuming spherical symmetry, let us again use m as the independent. vari- 
able, which is the mass enclosed within shell radius r. We define a new variable 
w, which is the volume of the sphere contained within r. We represent the full 
differential with respect to mass m as an overhead dot-symbol (), so that: 


ya ey på (59) 


_ dw 1 
Ww = — = = — 
dm p 

dw dv 


"dn dn 


The total energy Ftotal of the planetary body is a combination of its total 
internal energy Finternal,; Which include both the thermal energy and the elastic 
energy stored due to reversible-adiabatic (isentropic) compression, and its own 
gravitational energy Fsrav. We can write it in integral form as: 


Etotal =: Einternal gg Egrav (60) 
m=M, surface G.m 
= u— -dm 
m=0, center r 


Now, let’s introduce the Lagrangian L = L(m;w,w), with natural variable 


7099 


m and dependent variables (w,w), separated by the ”;” sign, as: 
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Then, the integral form of total planet energy Ftotal becomes a functional 
(function of function), which varies depending on the functional form w(m), and 
its first-order derivative w(m): 


m=M, surface 
Etotal = / —L(m;w,w)-dm (62) 
m=0, center 
m=M, surface G.m 
= u(w) — —— | : dm 
m=0, center hun sf 3 1 
Part I 3 -W3 
An 
Part IT 


Then, the problem of solving for planetary interior structure is reduced to 
finding the appropriate functional form of w(m) which minimizes the total en- 
ergy Frota] of the planet, provided that the specific entropy s = so within it 
remains constant. 

Part I in Eq. is the EOS of material. It only depends on w. It comes 
from the specific internal energy u = u(v, so) dependence on specific volume 
v = w and specific entropy s = so. So u is a uni-variate function: u = u(v). 
We then write out the derivative of u with respect to v and represent with the 
prime-symbol (*): 


u =u (v) =u (w) = — =-P (63) 
dP dP dm dP 1 


dv dm dv dm Ù 


S 
Il 
S 
RA 
e 
ma 
l 
S 
P 
& 
l 
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From another perspective, Part I (EOS) is a short-range interaction so it 
only depends on w, while Part II (gravitational interaction) is long-range so it 
only depends on w. 

We could then apply the one-dimensional Euler-Lagrange equation to Eq. 
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d OL OL 
Å Da ) - des en 


MY 
act only on Part I act only on Part II 
Recall that: 
OL Ou(w) 
— = — = P 
Ow Ow (63) 
Then, this gives rise to the following second-order differential equation: 
dP 4 me 
== ss 2s 
a u” (w) -ò= (=) G 473 (66) 


The negative sign (—) suggests the Pressure P is decreasing from center 
outward, which is what we expected. After removing the negative sign (—) in 
Eq. we can re-cast it into the following second-order differential equation: 


1/3 
dw dw AT 
4 ; = ~Gem-w4/3 
E (=) dm? (3) pe er) 


There are several important features of Eq. 


2 . . . . 
e $% > 0, since u” > 0. That is, w(m) is a concave in m. 
dw 1 
e — = . 
dm m=0 Pcenter 
d?w 


2 


am? |a 7 0. Close to the center, w(m) is approximately linear in m. 


Eq. [67] can solved, to obtain w(m), w(m), and &(m) as a explicit functions 
in m, by supplying the appropriate Internal-Energy-EOS: 


u = ulv, so) (68) 
and by supplying the appropriate boundary conditions: 


w(0) = 0, (69a) 
1 
“ = (69b) 
T m=M Psurface 


We can also write Eq. [66] in an integral form: 


1/3 
m | Ar G- m -dm 
P(m) = Poenter | (3 um) (70) 


Eq. [70|can also be understood as an Integro-differential Equation: 
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dw NT ™ f Ag ue G -m - dm' 
P (=) s%0 | = Poenter I 4/3 (71) 
dm o (81) ` fom) 
If the appropriate Pressure-Density-EOS is supplied as: 


Again, the goal is to solve for w(m) that satisfy Eq. [70] and the EOS Eq. 


2.2 Virial Theorem 


Recall that Euler-Lagrange Equation gives: 


d OL OL 
dm ( % J5 æ ex 
ʻa 
act only on Part I act only on Part II 


Multiply Eq. [73] with w on both sides, integrate from an arbitrary starting 
mı to an arbitrary ending məz, and substitute the independent variable in the 
integrant as m’, we have: 


m=m d (ƏL m=m aL 
. —|- dm! = -——. dm’ 74 
f, = dm’ (5) de f, j Ow “ ( ) 
The RHS of Eq. [4| can be written as: 
m'=mə aL m'=m2 OL 
adm = - dm! 75 
hi Ow å J. Olnw we ( ) 


The LHS of Eq. [74]can be integrated by parts as: 


m'=m2 aL acy |5 pm'=m aL 
Lev aia (5) = (5) m'=m E f. pa (5) (76) 
SS 


Part A Part B 


According to Eq. Part A of Eq. [76] is just the product of Pressure and 
Volume evaluated at the upper and lower limits: 


o(a) 


Part A 
Part B of Eq. [/6|can be re-written as: 


m'=m2 


= P(m2) - w(m2) — P(m1) : w(mı) (77) 


m'=m1ı 
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m/=m1q 


m'=m2 OL E m'=m2 ; dw OL m'=m2 å OL -{"™ OL 
f, g (3) 7 - NE ae Lo ve aa Ane 
= 
Part B 


(78) 
Therefore, combining all of these (Eq. [74 Eq. [75] Eq. [76] Eq. [77] and Eq. [78) 


into a single equation, we have: 


<S“ <S“ 


act only onw act only on w 
(79) 

Eq. [79] is the general form of virial theorem, applicable to any mı and mg. 
It is convenient when applied to a single layer that has a uniform composition 
within a multi-layered planet. At both the upper and lower boundary of this 
layer with adjacent layers, Pressure is continuous and Volume is continuous, and 
thus, their product P(m) : w(m) must be continuous. 

One special case is to take the integral from center m’ = 0 to surface m’ = M. 
P(m) : w(m) vanishes at the center since Volume is zero there. 
P(m) : w(m) vanishes at the surface since Pressure is zero there. 


Then we have: 


m'=M Ə Ə 
= J l - L(m';w,w) + dm" (80) 
m!=0 Olnw Olnw 
Yn RA 


act only onw act only on w 


If, however, there is certain pressure exerted onto the planet surface, due to 
for example the presence of an envelope, or an external field, then, 


P 7 7 i ð 4: 0 L(m';w Ww) dm! (81) 
surface * Vsurface = Er Olnw Olnw o? 
Sy ——" 


macroscopic microscopic 


As denoted above, part of the operator gives the macroscopic interaction, 
that is gravitation, and the other part of the operator gives the microscopic 
interaction, that is elasticity. 

The virial theorem can be viewed as a global relationship between the elastic 
energy and the gravitational energy of the bound object—planet. 

The form of Eq. invokes that perhaps we should re-cast the boundary 
conditions in a more symmetrical way as follows: 
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P(0)-w(0) =w- or Ma 0, (82a) 
P(M)-w(M)=w- a _ =0. (82b) 


2.2.1 Relation to Similarity Transformation 


This result (Eq. can also be viewed from the variational principle itself, by 
considering a small variation of the total action S about the equilibrium: 


m'=M 
0 = ôS = (Z iw GE 80] am (83) 
m'=0 ðw Ow 
Let’s consider a particular kind of variation: the similarity variation as 
dw = a-w. Here we take a as a small number (and a constant), then dw = a- w. 
It satisfies the first boundary condition of w(0) = 0 at m = 0 automatically. 
But it seems to violate the second boundary condition of w(M) = 1/ Psurface 
at m = M with this proportional variation or similarity transformation. 
However, noticing that the pressure P = oh is zero at the surface, so the 
effect of this variation vanishes at the surface (m = M) also. Therefore, by 
adopting this particular choice of dw, the same conclusion is reached: 


m'=M m'=M 
0=6S= (Hawt Eao) am =a f ( Z j 


m/=0 Ow Ow n/=0 ðlnw Olnw 
(84) 
From this perspective, virial theorem can be viewed as a special case or 
a direct consequence of the variational principle (stationary action principle) 
itself. 
The same argument of similarity transformation has been invoked by R.P. 
Feynman to derive the virial theorem for the generalized Thomas-Fermi model 
of atoms [40]. 


2.3 Normal Modes 


Now, if we consider the adiabatic radial motions (or pulsations) of a planetary 
body, we could also formulate it in the variational principle. This time, we need 
to include kinetic energy of motion in the Lagrangian £. 

Here, for convenience, we denote partial-derivative relative to mass oon as 
subscript m, and partial-derivative relative to time & as subscript +. 

Now, the volume w not only depends on mass m, but also depends on time 
t, to account for motions (in the radial direction): 


v=o mt (85) 


now two independent (natural) variables 
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Correspondingly, the Lagrangian £ becomes: 


pr 


now two independent (natural) variables 


c=e( m,t mer) (86) 


The full action functional S that we want to extremize with respect to 


w(m, t) is: 
Sw] = [J ehrt) -dm - dt (87) 


More specifically, in terms of the lower and upper limits of the integral: 


t=tp m=M 
Si] = if J £m, te un, -dm- dt (88) 
t 


Et, m=0 


where ta and ty are arbitrary starting and ending time moments. 
The variational principle is understood as: 


6S[w] = 0 (89) 


The specific kinetic energy 7 (in the radial direction) anywhere inside the 
planet, that is, energy of radial motion per unit mass: 


1 Or å 
(2) m 


Recall the definition of w earlier, we have: 
1/3 
3 
p= GE) nyt ls (91) 
An 
so, J can be expressed in terms of w, and its partial-derivative with respect 
to time t, that is, w; = Da 


9 2 
Fl (2) (92) 


le 


N| = 


The overall Lagrangian £ is classically defined as the difference between the 
kinetic energy term 7 and the potential energy term U (that is why initially 
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we keep a negative sign (—) in front of U, in order to be consistent with the 
definition here): 
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Now, the Euler-Lagrange Equation becomes: 


OL = aa) =)= (94) 


Ow Ot\ Our, Om \ Owm 


Eq. will eventually lead to a wave equation, where we can solve for its 
characteristic frequencies (eigenvalues) and characteristic solutions (eigenfunc- 
tions). 

Plugging in £ from Eq. [93| we have: 


2/3 —1/3 
oL 2? GE) a8 up (GE) Com-w 4/3 


ðw 27 Vår 3 Vår 
(95) 
ƏL 1. /3\? 
Du: = 9 . (+) 4/8 We 
ð (ƏL 4 FGVE oe 
(aa) = = (z) etg (GE) wn 
OL Ou(Wm) 


= = 1 = 
Owm T Own u (Wm) ~~ P(Wm) 


Om \ Owm m Owm 


Since the first partial-derivative of w relative to m is the reciprocal local 
density p(m) at m: 


1 
Wm = = = (96) 
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Then, we can see how u” is related with the material EOS properties, in 
particular, the isentropic bulk modulus K, and the isentropic sound speed cs: 


Plwm) _ OP(p) 3 P) APL) 
1 m 2 
m) = = — . = . = -Ks 
u (Wm) EYR Tee Ea (p) 
2 

K; Ks c 

p p Win 
Here, Ks = (2), is the local isentropic bulk modulus that we have 


defined earlier. cs = V= is the local isentropic sound speed in a fluid where 


shear modulus is zero (a fluid does not sustain shear stress). In general, both 
K, and cs are functions of density and can be calculated for a given EOS. 
Furthermore, u” is related with the second adiabatic index T: 


_ OP(wm) ə P (fOmPyV| — 
ua) Owm E p Olnp pee (28) 


In order to do that, let’s now apply the perturbation theory, and consider only 
small motions deviating from the stationary state. Suppose we can decompose 
w into two parts: the time-independent stationary part wo corresponding to the 
un-perturbed state (0-state), and the time-dependent perturbation dw which 
depends on both m and t. 

Furthermore, we apply separation of variables to separate the perturbation 
dw into a time-independent amplitude v(m), and a time-dependent oscillatory 
part e7: 


w(m,t) = wo + dw 
= wo(m) + dw(m, t) 


= wo(m) + (rm) - exp (ict) 


= wo(rn) (14 ct): exp (iot) ) (99) 


Generally speaking, both ¢ and o can be complex numbers, in order to 
accommodate phase shift in superposing different normal modes. 

a is the (complex) angular frequency of the oscillation. 

¢ is the (complex) fractional amplitude of the oscillation. It is dimensionless 
and |¢| < 1 for small perturbation. It vanishes at the center ¢(0) = 0 due to 
spherical symmetry. 

On the other hand, % has the dimension of Volume. 

Therefore, we can write out all the partial-derivatives of w in m and t as: 


wi = wo(m) : ¢(m) - (io) - exp (ict) (100) 
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2. exp (ict) 


(m): (io 
(m) : (—0°) - exp (iat) (101) 


IN IS 


Wm = 00,m T | Wo, m * ¢ + wo - 3 < exp (ict) 


| 


Wo: c) - exp (ict) (102) 


Wmm = W0,mm T (0m . ¢ +2- WO,m * Gm + Wo ° Som “exp (iot) 


= Wo,mm + (0 . c) - exp (ict) (103) 


Since |¢| < 1, we can neglect any term of the order Ç? or higher, thus, we 
don’t need to be worried about the w? term. Moreover, we know that wo(m) 
satisfies the stationary-state Euler-Lagrange Equation. 

By only keeping the term up to linear in ¢, we have (from Eq. [95): 


aos ea Gomer dett) ao 
slam) see 
2 (ZL) vr (dom + (20: Due): (sam (20° une) 
~- (e (wom) +u” (wom) (w0: Dan 5) . (sam + (wo: gn” se) 
= = (u (enm) mm + [u (nm) or Ga +" rm) (o0 Citam] =€) 
~ = (u (vam) sinom + [u (nm) 00:n], 1e) 


Now, let’s re-consider the stationary-state equality from Eq. 


Owm 


-1/3 
u” (wo,m(m)) - wo,mm(m) = É . ( 3 ) -Gem|- (wo(m)) 43 (105) 


år 


It can also be written as: 
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1 —1/3 —4/3 
v= [3 (2) Gem) å (106) 
2 
ar / (22) (107) 
1 3 pe 4/3 dwo ? d?wo 
a EEN (2) J (32) A 


By supplying an appropriate EOS, such as in the form u = u(v) = u(1/p), or 


P = P(v) = P(1/p), or cs = cs(v) = cs(1/p), the complete solution of Eq. 
will give all of the followings as explicit functions in m: 


e ulm] = u((wom(m))) 
e um] = u'((wo,m(m))) = -P[m] 
eum] = u” ((wom(m))) = (eslm])?/05,m 


after removing the stationary-state equality, Eq. [94] becomes: 


3 \-1/3 . 
0 GE) im 


4 
9 ÅT 
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3 A = io 
+5 (E) AE 


+ [w (wom) i (wo å Om] er (109) 


m 


Removing the common factor of e"7*, that is, completely removing the time- 
dependence and turn the equation into an ordinary differential equation only in 
variable m, and re-arranging the terms, now we have: 


d Wo ` ae Ne 
e OOE 
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Eq. is a Sturm-Liouville form. 
Now if we divide both sides Eq.[110|with u”, and substitute in the stationary- 
state identity of Eq. 


d d å mm å mm 
os (e (wo 2) | (p2 42 EEE (111) 


u" dm dm 3 Wo AT 


or equivalently, if we factor out womm/wo term: 


1 d d(w - mm 4 2 
EEG O) y (as pE) eoo E 


Now, we realize that the ratio m/wo is just the unperturbed mean density 
(p) enclosed within m. The term 47G7 reminds us of the dynamical timescale, 
and thus, we can define a new frequency parameter « as: 


_ /4nG mo AG > G-m 
a=) 3 way 3 p= rå vs 


k = &(m) decreases from center outward because p or p decreases from center 
outward. « is in fact the angular frequency of a circular orbit at radius r around 
mass m. 

Alternatively, one can also express wo(m) in terms of k(m): 


wlm) = FË. ær) (114) 


Therefore, Eq. [66] can be written as: 


8/3 
u” * WO, mm = oi = Å å am) (115) 


and, 


ul’ -w irtiri 3 k(m)!4/3 
TI ~ Ier? e) (116) 


We can re-write Eq. |112] as: 


d d . Wo,mm a 
båfe ta) (A) eB) ear um 


We recognize that (wo - Ç) is just the (complex) amplitude w of the volume 
perturbation defined earlier: 
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1 d „ dy Wo, mm 1 o? 
| $ . . 4 — . = 11 
u" dm (2 dm wo 3 K K? y= 8) 


If we multiply both sides Eq. with (u”) to get back to the original 
equation, we have: 


d | „d u"-womm \ 1 g” 
; Fr FEN VE ti 
dm (2 w) = ( Wo 3 a k? ped vi 


which, in the language of k, becomes: 


d | „dẹ 1 k(m)!/3 o? 
: -{4+ ——~ | :v=0 120 
dm (2 dm i 167? \ (G- m)4/3 i k(m)? y ey 
The advantage of expressing everything in « is that « is positive finite and 
monotonically decreases (slowly) from center outward. u” = pPT = pK, = på 


is positive finite and monotonically decreases (slowly) from center outward. 
Consequently, we have from Eq. and Eq. 


1 d | „dy 1 (m) !4/3 o2 E 
u” dm ( w) T Toru" ` Cee Ve ey) 


has dimension of 1/Mass? 


We can define a new parameter À. To the leading order, Å ~ m2/3. Thus, we 


(m) | 44 22 2 


PIE) 


can write: Mm) = a(m)-m?/3, where a(m) = 47 vu". Or 1 


where a(m) is positive finite. Near the center of planet, where both u” and « are 
approximately constants, then a(m) const, which has the dimension of Mass: 


Mm) = 47 va" > (Gant) E ) (122) 


k(m) 1442 


k(m)? 


Then, Eq. can be re-written as: 


1 d | „dọ 1 - 
u” dm (: =| a A(m)? Ke (12a) 


One can think of Å as the (local) wavelength in the mass dimension. That is, 
after ~ A/2 in the m-coordinate, the volume perturbation dw changes sign from 
(—) or (+) or vice versa. Recall that Vu” = pc, where p is the local density, and 
cs is the local effective sound speed, which measures how fast a perturbation 
(either in density or pressure) signal propagates in a medium. Then, 
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À (G -m)?/3 1 
= = 2T pcs * : (124) 
Å n(m)? LF 


Remember that o can be a complex number, in particular, if it is imaginary, 
i.e., multiple of i, then it means exponential growth (multiple of —i) or decay 
(multiple of +i) in the system (see Eq. [99h. That suggests the possibility of 
emergence of structure in the system. 

However, if ø is purely imaginary, |o| cannot exceed 2k, in order to make 


4/4 aay well-defined. This can be interpreted as that the gravitational 


growth or collapse cannot happen faster than the dynamical timescale, as we 
expect. 

One special situation is when ø = 0, then the perturbation is stable and 
neither grows or decays nor oscillates with time. When this happens: 


(125) 


This may be interpreted as the characteristic mass units that emerge out of 
the total system mass m of a gravitationally-bound system with signal speed 
Cs, if the whole system is only slowly evolving and structures emerge slowly, 
compared to the dynamical timescale of (1/«). It may help explain the formation 
of ring-like structures in proto-planetary disk, as those structures are stable over 
at least multiple orbits. Our derivation is general. 

If things are sticky on the small-scale, then it may help form structures on 
the large-scale, especially true in the molecular clouds when the initial density 
is small. Stickiness is manifested by the molecular forces of cohesion and can be 
characterized by a negative internal pressure term. This cohesion would help 
form structures. Van der Waal force, though not very strong, acting over a large 
ensemble, may become significant. That Vu” in Å may as well be V—u”, and 
then the sound speed c, represents how fast the stickiness propagates. :) 

Generally, we should think of u” as some kind of coupling at local scale 
that either resist or assist collapse. The critical molar density, for the stickiness 
(coupling) to take effect, as we will show later, is ~ 0.01 mol/cc (or ~ 0.001 
mol/cc at lower temperature). Thus, structures not only emerge from top down, 
but may also emerge from bottom up. 


2.4 Derivation of Wave Equation 


If we multiply both sides Eq. with (u”)?, we have: 


„d „ aw u”? 
i (: Je ;-v=0 (126) 


dm dm 


Now we can define a new coordinate x to replace m: 


36 


dx = dm/u", 


fo a (127) 
Ju" fr J K; 


Then, Eq. |126| transforms into a differential equation in x = f dm/u": 


dy 1 
da? * Dune = oe) 

Eq. [128] resembles the time-independent Schrödinger Equation. (A/u”) can 
be thought of as the new (local) wavelength measured in the new coordinate zx. 
Generally speaking, (A/u”) is not constant and varies with x. 

The new coordinate x actually has its physical meaning: dæ is the volume 
element dw modulated by the local bulk modulus K,, in the un-perturbed state 
(0-state). It is the natural coordinate viewed by a wave propagating in the 
interior of a planetary body. In future, we will emphasize the importance of 
natural coordinate, and when using it, greatly simplifies a physical problem. 
In fact, the least action principle with its action S and Lagrangian £L can be 
expressed in whichever coordinate more convenient for calculation. 

x has the dimension of Volume/Pressure. It ranges from 0 to £m 


m=M 
tm = i on (129) 


Then, the whole problem is reduced to solving this wave equation (Eq. |123 
or Eq. for its characteristic frequencies å (eigenvalues) and corresponding 
characteristic solutions w(m) or y(x) (eigenfunctions). 

o is quantized, that is, ø can only take discrete values, because of the bound- 
ary conditions. 

For Eq. [123] the boundary conditions are: 


w(m = 0) = 0, (130a) 
as i (130b) 
dm 

m=M 


For Eq.[128] the boundary conditions are: 


p(x =0)=0, (131a) 
dv 7 
dr i =0. (131b) 


The first boundary condition (Eq. [130a] or Eq. [131a) occurs because of spher- 
ical symmetry, that the amplitude ~ must vanish at the center of the planet. 

The second boundary condition (Eq. [130b] or Eq. [131b) occurs because the 
gradient of pressure disturbance which drives the differential volume perturba- 
tion vanishes at the surface of the planet. Although not immediately clear, this 
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condition is equivalent to requiring that all interior disturbances be reflected at 
the surface (as itself moves) back into the interior; that is, no pulsation energy 
is lost from the interior into space because all is reflected back inward from the 
surface [46]. See also in [58] [81] [32]. 

Generally speaking, Equations of the form of Eq.[118] [H8]and Eq.[128Jare difficult 
to solve analytically, because the coefficient in front of y is not a constant but 
varies with m or x, which means the wavelength varies with location. 

Even the following simple example requires the introduction of a special 
function called Airy function (introduced by the British astronomer G. B. Airy): 


2 
= +x.v=0 (132) 
da? 
Historically, many techniques to solve such equations have been developed, 
such as the ladder operators or the WKB-approximation. 
Indeed, if |A/ u"| < zm, that is, the wavelength is small compared to the 
overall dimension, then ~ can be approximated as: 


va) EE (133) 
Or, 
add i f dm/u" ve —i f dm/u" 
Vix) A-e PA ur) +B VE (134) 


This is the essence of the WKB-approximation. 


2.5 Legendre Transformation 


This product: (w: P), is important. Using it, one can define another equivalent 
action S and equivalent Lagrangian Å in terms of using Pressure P as the 
dependent variable, instead of volume w. The indepedent variable is still mass 
m, as follows: 


L+L= 7 (Ww P)=w-P+a-P (135) 
in other words, 
~ d 
fe Lø Pl (136) 
dm 
so that, 
Stø) - fea KE 
w| = m, dm Ww ’ 
» at ox z 
sip)= fe din; Í spl = 
dm 


38 


Here È is the Legendre transformation of £, where the variables become P 
and P as: £L = L(m; P, P) versus £ = L(m;w,w). 


Moreover, 
dL+dL=w-dP+P-dw+w-dP+P-dw (137) 
Recall that: 
L(m;w,w) = —( u(w) = fo) m ) 
A —— 
close-range interaction due to EOS  far-range interaction due to gravity 
(138) 


Here f(w) encodes the geometric shape of the system. In particular, for 


spherical symmetry, 


Fu) =G/r=C- 2) T (139) 


2.6 Hamiltonian 
Hamiltonian H is also related with Lagrangian £, but with the transformation 
carried only half-way through. 


(140) 


where, h(P) is exactly the definition of specific enthalpy, assuming (1) the 
process occurs along an isentrope, and (2) no phase transition occurs. 
The Symplectic Equations can now be written as: 


P= ar = f'(w)-m, (141a) 
os “ =h'(P), (141b) 
OH aL 

am ag = Tw) (141c) 


Eq. is in fact the hydrostatic equilibrium equation. Eq. |141b}is in fact 
the EOS. Eq. is the mass-evolution of Hamiltonian—Symplectic Integral! 
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2.7 Hamilton’s Principle Function 


From the perspective of Hamilton-Jacobi Theory, the above formulation can 
be further distilled and summarized by introducing the Hamilton's Principle 
Function. The Hamilton’s Principal Function is defined as the function of the 
upper limit of the action integral taken along the minimal action trajectory of 
the system: 


(m,w) 
S(m,w) := / L-dm (142) 


This is to view to total action S as a function of m and w on equal footing, 
with the upper limit of the integral variable. In this form, S does not explicitly 
depend on either w or P. Calculating the variation of with respect to the 
endpoint (m,w) gives: 


OS 
AG = P (143a) 
OS 
Am = TH (143b) 


Then, the Hamilton-Jacobi Equation (HJE) is a single, first-order partial- 
differential equation based on the Hamilton’s Principle Function S(m, w), which 
is equivalent of an integral minimization problem such as the Hamiltons Prin- 
ciple, and can be used to determine the Geodesics. 


OS(m,w) _ _ dS(m,w) 
ae (mie, a ee 
Plugging the particular form of Hamiltonian H (Eq. considered here: 
- OS(M,w)N ,(OS(m,w) 
Then, we have: 
å OS(m, w) ai Mg (146) 
Ow om 


Ultimately, we want to solve for the Hamilton’s Principal Function S(m, w) 
using Eq. [146]given that h (EOS) and f (geometry and symmetry) are provided. 

If we express S in radius r instead of volume w, then, HJE can be written 
as: 


=0 (147) 


Amr? Or r om 


il 1 -o Gem , 3S(m,r) 
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2.7.1 An Example 
For EOS of: 
P=K. pP (148) 
Then, the specific enthalpy h(P), by definition, is 
h(P) =2-VK.VP (149) 


Let’s define a new parameter: 


27G 
k= V= (150) 


k has the dimension of 1/distance. 
Then, the Hamilton’s Principal Function S(m,r) for this EOS is: 


ls ( om) l (am =) vey) 


According to Kq.|143a| the Pressure is: 


_ OS(m,w) 
7 Ow 
1 dS(m,r) 


~ fnr? Or 


E zo l (Å ou) (art) 
= yan (£ = a?) l (i r) oe sin (k - 5 | 


| 7 


P 


Then, 


h=2-VK-VP 
-2 Ore [Em (ene) 


22. (a) å 


Then, 
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(1-9) Ears) 
k 


=G-m- | ———=— 154 
ae (Ga) es) 
On the other hand, 


OS(m,r) k 
3m SEE (arta) vee 


Plugging back into Eq. one can see it is satisfied. 

The essence (spirit) behind this method is Power Series, as Cauchy used to 
solve Kepler Equation: E — e - sin E = M [66]. In a sense, tan(k-r) —k-r can 
be understood as a power series expansion of radius r around the origin: 


17 
315 


(kr) 4 n (k-r)? +... (156) 


2 
2 ihe 
+ 7p (her) + 2835 


1 
(kr) 
(kryt+g 


tan(k-r) —(k-r) = 5 


Or alternatively, in terms of Laurent Series: 


1 3 6 (k-r) Ak-r) 
tan(k-r)—(k-r) (kr)? 5(k-r) 175 7875 
37(k-r)> 118(k-r)"  5506(k- r)? 
3031875 197071875 186232921875 


+ 


(157) 


Therefore, in this particular case, S(m,r) can be expressed as: 


G- m? 3 6 (k-r) 2(k-r)8 
S(m,r) = (- =) “ke (a 5(k- 1) 175 7875 


37(k+r)? _118(k:r)"  5506(K re | (158) 
3031875 197071875 186232921875 |” 


In particular, when (k - r) = 1/2, then, S(m,r) = 0. 
In particular, when (k - r) = 7, then, S(m,r) = En Å, 
Since only partial derivatives of S are involved, there could be an arbitrary 


additive constant to S. 
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3 Hydrogen EOS 


Hydrogen EOS is important as it is relevant to the interior of many gaseous 
planets. It has been subject to intensive experimental investigation [67]. 


3.1 Critical Point (C.P) of Liquid Hydrogen 
The Liquid Hydrogen (H2) has its critical point (c) as follows [10]: 


pc ~ 0.03012 g/cc 
Te = 33.20 K = —239.95°C 
P. =1.3-10° Pa = 13. bar 
ne ~ 0.03012/2 ~ 0.015 mol/cc 
(159) 
Here ne is the molar density expressed in (mol/cc) of Hə at its critical point. 
Te of Hp is very low and hence the hydrogen in planet interior should be consid- 
ered super-critical molecular fluid when it is above the molar density of 0.015 


mol/cc. At even higher density or temperature, it then transitions into metallic 
fluid. 


3.2 Tripe Point (t1) of Liquid Hydrogen 


The Liquid Hydrogen (H2) has its triple point (t1) co-existing with solid-H> and 
gas-H» as follows [10] [6] [51]: 


Pti,Liquia © 0.07691 g/cc 
pt1,Solia © 0.086 g/cc 
Ta = 13.96 K = —259.2°C 
Pi = 7.7- 10° Pa = 0.077 bar 
Nt1 Liquid © 0.07691/2 ~ 0.0383 mol/ce 
N#1,Solia © 0.086/2 ~ 0.043 mol/ec 


(160) 


Here nt ,Liquia is the molar density expressed in (mol/cc) of H2-liquid at its 
triple point. Here n¢1,golia is the molar density expressed in (mol/cc) of Ha-soild 
at its triple point. The solid-Hz was obtained by James Dewar in (1899) and is 
one of the lowest-density solids [35]. 

See also in Fig. [6] It is interesting that: 


mi/nev 3 (161) 


We will re-visit this point in a later discussion and generalize it to other 
molecules. 
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3.3 Mazevet 2019 EOS 


We apply this technique to Mazevet 2019 Hydrogen-EOS [29]. First, we visual- 
ize this EOS in 3D the logarithmic relationship between pressure, density and 
temperature (Fig[8). 

This 3D EOS can be projected onto various cardinal directions, to show 
the T-p (Temperature-Density) relation (Fig[p) or the P-p (Pressure-Density) 
relation (Fig[10) along these isentropes. On the other hand, the projection onto 
P-T (Pressure-Temperature) plane is often encountered in illustrating a phase 
diagram in the literature. 

With this understanding of EOS, one can then convert these isentropes into 
corresponding mass-radius curves (Fig.[11), by solving the Differential Equations 
mentioned in the previous section (Mathematica code is provided in a separate 
file). These mass-radius curves are for envelope only, not including the atmo- 
sphere. The calculation is truncated as density drops below 0.01 gram/cubic 
centimeters. There are two reasons for this truncation: 


e The atmosphere temperature profile is likely more isothermal due to the 
stellar irradiation, as opposed to isentropic that is driven by interior con- 
vection. In other words, the specific entropy of the envelope is nearly con- 
stant, and is determined by the outward energy transport in the planet 
interior, while that of atmosphere is influenced strongly by the host star. 


e The atmospheric scale height of the planet is relatively small, on the order 
of 0.1 Earth radius, for H> atm of 1000K. If considering 10 scale-heights, 
it would be a correction of 1 Earth radius. We could easily add that later. 
The atmosphere can sometimes be treated with modified ideal gas EOS. 


Furthermore, one can add contours to illustrate the underlying physical 
meanings of these mass-radius curves. For example, one can add two sets of 
contours, one set for specific entropy S (eV/1000K/atom), and the other set for 
central density of planet in po (g/cc). When a planet cools but not losing mass, 
it moves vertically (because mass is conserved) on the mass-radius diagram. 
From this point, one can get a sense of how its interior contracts (with central 
density and bulk density increasing) over cooling. (See Fig. 


3.4 Becker 2014 EOS 


Now, let’s perform a detailed comparison between Mazevet 2019 EOS with that 
of Becker 2014 EOS from the German Rostock group ([I7], referred to as H- 
REOS.3, where 3 stands for version 3). 

For Becker 2014 EOS, again we construct the temperature-density relation 
(Fig[9) or the density-pressure relation (Fig[10) along eight isentropes of differ- 
ent specific entropies. 

Carrying on the same procedure of solving differential equations, we obtain 
the mass-radius curves corresponding to Becker 2014 Hydrogen EOS (Fig. [15). 

Then, we make a comparison between the Mazevet 2019 EOS and Becker 


2014 EOS (Fig. [16]. 
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Figure 8: Hydrogen EOS of Mazevet 2019 with isentropes visualized in the 
three dimensions of lgT, lgp, and lgP. The eight color curves corresponding 
to isentropes of specific entropy S (eV/1000Kelvin/atom) = 0.3, 0.4, 0.5, 0.6, 
0.7, 0.8, 0.9, 1.0. The bulge and the oscillation behaviors of the isentropes 
with lower specific entropies as a result of this bulge, are due to the Plasma 


OlnT 
Olnp 


Phase Transition (PPT). Recall that along an isentrope, y = ( ) , and 


sS 


r= (317 ) . Therefore, y and I are just the tangential slopes of these 3-D 
sS 


isentropes projected onto the T-p and P-p plane correspondingly. 
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Figure 9: lgT-lgp relation of isentropes in the Hydrogen EOS of Mazevet 2019. 


The tangential slope in this plot is y = (382) . We call it a y-plot. 
S 


3.5 Discussion 


The major difference between Mazevet 2019 EOS and Becker 2014 EOS are their 
treatment above 10° Kelvin, resulting in the biggest difference among the mass- 
radius curves of the highest entropy. Specifically, Becker 2014 EOS does not 
have a jump in thermodynamic properties such as entropy at 10° Kelvin, or at 
density of 10 g/cc, while Mazevet 2019 EOS does. These jumps in Mazevet 2019 
EOS are artificial, because the authors stitch together different EOSs calculated 
from different methods, applicable to different lgT-lgp regimes (From Figure 1 
of [29], 10° Kelvin is boundary between SCvH and H-EOS, and 10 g/cc is the 
boundary between QMD and CP98.). 

This is better visualized in 3D the relationship between specific entropy, 
logarithmic density and logarithmic temperature (Fig[17). 

The three hottest isentropes (the three right-most curves) go across a jump 
at 10° Kelvin, into the so called ” Dragon” regime in a generic EOS diagram 
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Figure 10: lgp-lgP relation of isentropes in the Hydrogen EOS of Mazevet 2019. 


The tangential slope in this plot is I = (Re ) . We call it a I -plot. 


(Fig[18). That is crossing the boundary of Fermi-energy equals Thermal energy 
of electrons. This ” Dragon” regime is complicated because it is both partially 
degenerate and partially ionized, so the EOS there is not very well understood. 
It would require fully quantum numerical calculations such as Path-Integral 
Monte Carlo (PIMC) or Path-Integral Molecular Dynamics (PIMD). 

According to Dave Stevenson [79], the interiors of many gas giants are simple. 
To a very good first-order approximation, they can be characterized by ONLY 
two parameters: 


e central density po 


e either the planet radius r, or its interior specific entropy s (constant 
throughout its interior) 


The reason can be seen in the mass-radius curve grid itself: If one looks at the 
mass-radius diagram corresponding to Mazevet 2019 Hydrogen EOS (Fig. 
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Figure 11: Mass-radius curves of Hydrogen EOS of Mazevet 2019: the same 
color schemes of isentropes apply, but do not correspond to the temperature 
color scale on the upper-left. The temperature color scale is only for exoplanets 
plotted. 


or Becker 2014 Hydrogen EOS (Fig. 15), the constant central density contours 
are almost parallel to the constant bulk density contours (the gray dashed lines 
in the background marked by 0.3, 1, 3, 10 g/cc). This implies that the interior 
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Figure 12: Modified Mass-radius curves of Hydrogen EOS of Mazevet 2019 


structures of these planets are homologous to one another. The density every- 
where inside a hydrogen planet almost scales linearly with its central density. 
In most of the parameter space, the isentropes in the pressure-density (P-p) 
parameter-space can be very well approximated by the parabolic power law 


relation of P ~ på (See Fig. [19). 


This can be better understood by the following two examples: Fig. 
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Figure 13: lgT-lgp relation of isentropes in the Hydrogen EOS of Becker 2014. 


The tangential slope in this plot is y = oe . This is the y-plot. 
S 


Fig.|20b| The oscillation behavior is general for any cubic or higher-order poly- 
nomial EOS at phase transition, and can be solved by Maxwell construction. 
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Figure 14: lgp-lgP relation of isentropes in the Hydrogen EOS of Becker 2014. 
“ine | . This is the I-plot. 


The tangential slope in this plot is I = EW 
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Figure 15: Mass-radius curves of the Hydrogen EOS of Becker 2014 
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Figure 16: Mass-radius curves comparison: Mazevet 2019 EOS (dashed) versus 
Becker 2014 EOS. 
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Figure 17: Hydrogen EOS of Mazevet 2019 with isentropes viewed in three 
dimensions of lgT, lgp, and S (specific entropy in unit of eV/1000Kelvin/atom). 
The eight color curves corresponding to isentropes of S = 0.3, 0.4, 0.5, 0.6, 0.7, 
0.8, 0.9, 1.0 . The bulge and the oscillation behaviors of the isentropes with 
lower specific entropies as a result of this bulge, are due to the Plasma Phase 
Transition (PPT). 
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Figure 18: Generic EOS diagram. Copyright 2011 by Alberto Douce [36] 
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Figure 19: black dashed lines correspond to P ~ p? in lg-lg plot. The tangential 


OlnP 


slope in this plot is I = alap 


This is again the I'-plot of Hydrogen. 


This plot suggests that I æ 2 for fluid metallic hydrogen over many orders-of- 
magnitude in parameter-space. 
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(a) an example of a hot super-Jupiter 
(1800 Mg): note that its interior den- 
sity is almost linear in depth z = (R-r) 
except at the center, and its interior 
pressure profile is quadratic in depth, 
and its temperature profile is almost 
linear in depth as well. This is all due 
to P ~ på being approximately true for 
metallic fluid hydrogen, and can thus 
be generalized to all other hot jupiters. 
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(b) Jupiter (317 Mø): Its tempera- 
ture profile shows an oscillation behav- 
ior similar to a van der Waals EOS 
isotherm crossing a phase boundary. 
Here it is due to the metallization of 
Hydrogen, also called the PPT-Plasma 
Phase Transition. In reality, the actual 
temperature profile should be replaced 
by something similar to a Maxwell con- 
struction. 


P (GPa) 


T (Kelvin) 


4 H:O EOS 


Recent observations have suggested the possible existence of water-rich exoplan- 
ets [87]. Thus, the study of H2O EOS and its application to the interior 
structure modelling of these exoplanets is essential. 


4.1 Fluid H,O: 
(Gas+Liquid+ Supercritical) 


The International Association for the Properties of Water and Steam (IAPWS) 
constantly integrates new experimental data and constantly publishes the EOS 

of H20 in Liquid and Vapor Phases, in tabular forms and also in analytical 
formulae through interpolation. 


4.1.1 Critical Point C.P. of H20: 
Principle of Corresponding States 


It should be noted that J. van der Waals’ idea of corresponding states was 
anticipated by D. Mendeleev (1870), who pointed to the usefulness of comparing 
volumes not at the boiling points, but at temperatures when the cohesion of the 
liquids is close to zero . This idea is based on the principle of corresponding 
states [53]. According to this principle, the same or similar common functional 


relationship: 
p T. P 
£ EI N= 162 
(LE) (162) 


exists for similar substances. In other words, instead of expressing the prop- 
erties of a fluid (liquid phase + gas phase) through p, T, and P, we now use 


dimensionless units for this purpose — the reduced parameters: 6 = å fet 


T= aa then we can obtain a reduced “universal” EOS f(6,7,7) =0. The goal 
is one arrow for a thousand birds. 

Many important corollaries follow from this principle. For example, if any 
two fluids, say H20 and NH3, have identical values of m and 7, then they 
should occupy an approximately identical reduced density 6. However, an EOS 
characterized by only three constants (pe,Te,Pe) cannot be infinitely accurate. 
Therefore, both the principle of corresponding states and the conclusions follow- 
ing from it are only approximate. The errors in calculating various properties 
typically do not exceed 10%, especially for homologous substances or substances 
with close boiling points, such as H20 and NH3. Fundamentally, it results from 
the similarities in the microscopic structures between H20 and NH3, that both 
are polar molecules with hydrogen bonds. Moreover, the bond strengths of the 
covalent bond O-H and N-H are similar, ~ 100 kcal/(mol-bond), or equivalently, 
~4.5 eV/bond. 

The analytical formulae of H2O are thus anchored at its Critical Point C.P.. 
Here underscript c to denote various parameters at critical point [60]: 
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Pc = 0.322 g/cc 

Te = 647.096 K = 373.946°C 

P, = 2.2064-10" Pa = 220.64 bar 
Ne = 0.322/18 = 0.018 mol/cc 


(163) 
Here ne is the molar density in (mol/cc) of water at its critical point. 
Reduced density 6 (dimensionless), 
= £ (164) 
Pe 
and inverse reduced temperature 7 (dimensionless), 
Te 
== 165 
r=3 (165) 


are used as independent variables, to express other thermodynamic variables 
and functions. The reason to use the inverse reduced temperature 7 is because 
generally speaking, the higher the temperature (the smaller the 7), the closer 
the substance behaves as ideal gas, and thus, the series expansion based on ideal 


gas EOS is better expressed in 7 rather than (+) 


4.1.2 Series Expansion based on Ideal-Gas EOS 


All other thermodynamic quantities can be derived from the the specific Helmholtz 
free energy f. (We follow our derivation here as in IAPWS) [13] [49] 
Fundamentally, the specific Helmholtz free energy f is: 


f=u-T:s (166) 
It can be cast into a dimensionless form as: 
b= F =Y, r) +6"(6,7) (167) 
R i T pl H 


w° is the ideal-gas part, w” is the residual part. 
w° and w" is each expressed as analytic formula in terms of å and 7, with 
coefficients determined by best-fit to experimental data. 


8 
wd, 7] = Ind+n{+n§-r+n§-Inr+ one -In[l—exp(—7?-7)] (168) 
i=4 


where n? and y? are dimensionless coefficients. 
Likewise, 
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7 51 
P" (6, 7] =Y ng t r" 4 Vi ng 68 ot exp (69) + 
i=1 1=8 


54 
JO må 6% -r exp (0 + (6 = €i)? — Bi (7 9) + 


i=52 
56 
Xni A-S (169) 
i=55 


with all kinds of coefficients introduced to fine tune the formula to better 
match the experimental data: 


A = 6? + Bi- [(6— 1)7|* 
O= (1-7) +A: (6-19)? 
w = exp (0; - (6 — 1)? — Di- (7 — 1)”) 


(170) 


A specific gas constant for this mass-based formulation is adopted: 


R = 0.46151805 kilojoules - kg™t - K~! (171) 


Much of the complexity is introduced to describe the behavior of H20 around 
its Critical Point (C.P.), where the pressure changes tremendously over a very 
small density and temperature range. Beyond C.P., the distinction between 
the liquid and vapor phases vanish, and water is supercritical, existing as small 
but liquid-like hydrogen-bonded clusters dispersed within a gas-like phase [31]. 

The pressure P can be calculated from the specific Helmholtz free energy f 
through fundamental thermodynamic relation: 


of 
P= >. (54) 172 
Er Å (172) 
and be expressed in dimensionless variables as: 
P(6,7) ay" 
= 1 . 1 


Likewise, all other thermodynamic functions and variables can be derived 
from the specific Helmholtz free energy f through appropriate partial-derivatives. 
Internal energy (specific) u: 


wear (Ge) 4G) e 


Entropy (specific) s: 
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PE) ow 


Enthalpy (specific) h: 


n(5,7) _ aur) (au) \ (avr 
ep ater (2) +(F)) ne (Gr). a) 
4.2 ‘Triple Points (T.P.) 


Triple Points (T.P.) are important anchoring points to make the p-T phase 
diagram. We use underscript (t) plus a number to denote different T.P.. 


4.2.1 Gas-Liquid-Ice Ih T.P. (t1) 


Furthermore, specific internal energy u and specific entropy s of saturated liquid 
at the Gas-Liquid-Ice Ih triple point have been set to zero to determine the con- 
stants of integration. The Gas-Liquid-Ice Ih triple point (t1) has the following 
properties [60]: 


pu = 0.99979 ~ 1 g/cc for liquid 
Ta = 273.16 K = 0.01°C 
Pi = 611.655 Pa ~ 0.006 bar 
my Z 1/18 ~ 0.0555 mol/ce for liquid 
(177) 


Since the entire formulation is based on Temperature T and density p as 
independent variables. We ought to construct a T-p Phase Diagram to visualize 
the behavior of water. 

The range of IAPWS formulation is the entire stable fluid (including liquid, 
vapor, and supercritical fluid beyond C.P.). It is experimentally verified for 
Temperature T in between 273 K and 1273 K, and Pressure P up to ~1 GPa. 
Tests have shown that this IAPWS analytical formulation of the H20 EOS can 
be safely extrapolated for density and enthalpy of undissociated H20 up to 
~5000 K and ~100 GPa at least. 

The isentropes (the adiabatic curves which most likely resemble the T-p pro- 
files in the interiors of H2O-layer in a planetary body) are calculated and shown 
here (Fig. pi). All the isentropes start from the Liquid-Vapor phase boundary 
(the solid green curve), except the one above C.P. that are calculated from 
the isochor of p=0.322 g/cc. The dashed part of the curves are extrapolation 
from the experimentally verified regime, all the way up to the Fluid-Solid phase 
boundary. Here, the Fluid-Solid phase boundary is the boundary of Fluid-Ice 
VII and Fluid-Superionic Ice. 
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4.2.2 Fluid-Ice VII-Ice XVIII T.P. (t2) 


The superionic ice is also called Ice XVIII (Ice-eighteen) [30]. It is now exper- 
imentally verified [63]. The triple point (t2) of Fluid-Ice VII-Superionic Ice is 
another important anchor point for plotting the phase diagram, with approxi- 
mate temperature of ~1000 K and pressure of ~50 GPa. 


Pia ~ 2.6 g/cc for fluid 

T2 ~ 1000 K 

Pia ~ 50 GPa ~ 5 - 10° bar 

ni2 © 2.6/18 = 0.144 mol/cc for liquid 


(178) 
4.2.3 Ice VI-Ice VII-Liquid T.P. (t3) 
Tig = 355 K 
Pig = 2.216 GPa = 2.216 - 10* bar 
pvr = 1.31 g/cc for Ice VI 
pvu = 1.567 g/cc for Ice VII 
PLiq = 1.35 g/cc for liquid 
ny ~ 1.31/18 ~ 0.0728 mol/cc for Ice VI 
ny & 1.567/18 ~ 0.0871 mol/cc for Ice VII 
NLiq © 1.35/18 ~ 0.075 mol/ce for liquid 
(179) 


where n is the molar density of each phase in (mol/cc) at this triple point. 


4.3 Solid: 
Various Forms of Ices 


Each form of ice requires its own EOS, because of the difference in crystal 
structure and compressibility among the different ice forms. There exist solid- 
solid phase transitions among the different ice forms. 


4.3.1 General EOS of Ices 


Ices can assume different crystal structures, depending on the temperature and 
density (or pressure) conditions. For a crystal, its microscopic structure dic- 
tates its macroscopic properties. The Ice Ih (Ice-one-hexagonal) is the common 
form of ice found in ambient conditions on Earth surface. Its low density is 
due to loose packing, which is a result of the orientation of hydrogen bond- 
ing. However, in the deep interior of planet, subject to tremendous amount of 
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Figure 21: Isentropes of Fluid-H20 in T-p space. C.P. stands for the Critical 


Point between Liquid-Gas-Supercritical Fluid. Dashed curves are extrapolations 
from experimentally validated regimes. The tangential slope in this plot is 


y= (322) . This is the y-plot of water. 


Olnp 
s 


pressure and heat, the ice crystals change its packing to become more symmet- 
rical and denser. One common form of such a high-pressure ice phase is Ice 
VII (Ice-seven). It is believed to exist in the interiors of icy satellites such as 
Europa and Ganymede. Recent theoretical and experimental works suggest the 
de-localization or mobilization of protons in ice crystals above a certain pressure 
threshold. The free protons (hydrogen ions) within the lattices of oxygen ions 
can flow and give rise to moderate amount of electric conductivity. This new 
phase is thus termed the Superionic ice phase [30]. Although the protons are 
mobile and fluidic, it is still considered a solid phase, because the oxygen ions 
are fixed in the crystal lattice. 
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104 


The pressure-EOS of any ice, and in general, for any solids, can be con- 
structed as follows: 


Protal = P-tastic (p) F P;hnermal(P, T) + Peiectiónie 0: T) (180) 


e Plastic is the pressure resulting from cold elastic compression, and is a 
function of density p only. 


e Pinermal IS the lattice vibrational contribution, and is a function of both 
density p and temperature T. 


e Plectronic comes from the heat capacity of electrons calculated from finite- 
temperature Fermi-Dirac distribution, and is a function of both material 
density p and temperature T. In this context, we ignore it for now. 


4.3.2 Born-Mie EOS 


A general form of Paastic called Born-Mie EOS is constructed from a modified 
power law, with one attractive term and one repulsive term: 


a = 
a as) 
Po Po 
where m and n are constants. 


Note that the pressure increment dP in Born-Mie EOS can be written as: 


3:K 
Pporn-Mic(P) = 2 : 


m— mn 


a OCO (182) 


Thus, the bulk modulus K at any density can be written as: 


m+3 n+3 
3 Po 3 Po 


(183) 
The bulk modulus increment dk in Born-Mie EOS can be written as: 


2 ar 2 2 
GER 
3 Po 3 Po Po 

Thus, the derivative of K with respect to P, can be written as: 

eat Gents 

K 
KE Me i o (185) 
dP m+3 (2) n+3 (2) 
3 Po 3 Po 
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3- K 
dPBorn-Mie(P) = 2 i 


m— mn 


i dP 38 Ko 
Born-Mie = dln (p/po) = vi 


_ 3+ Ko 


m— mn 


dk 


The asymptotic behavior of K'(= dK/dP) (dimensionless) is as follows: 


dk 3 3 
eZ) MY H+) (186a) 
P=P0 3 
dK 3 
p—> +00 3 


Therefore, K'(= dK/dP) is a decreasing function in p or P. In principle, 
both Kó and K% can be determined by experiments. 
In some cases, it is convenient to express m and n in terms of Kj and Kh: 


n=3-(Kj- Ki) — 3; (187a) 
m=3-K/, -3. (187b) 
Then, Born-Mie EOS can be re-written as: 
Ko p K p (K-K) 
P orn-Mie = . 
Botne) (2: Ko- Ko) (2) (2) Gi 


Thus, it is parametrized by three parameters (Ko, Kj, and K} ) to com- 
pletely determine its functional form. 

Sometimes, the number of required parameters can be reduced to two, if one 
takes K% = 2, and Eq.|188}can be re-written as [48]: 


anno = å (2) (i) | i 


Eq. only needs two parameters (Ko, K6). 


4.3.3 Birch-Murnagham 2nd-order (BM2) EOS 


For our purpose, a more specific analytic formula called the Birch-Murnagham 
2nd-order (BM2) EOS is adopted by setting m=4 and n=2 , to approximate 
Patastic; due to its compliance with experimental high-pressure data of mineral 
compression, and good asymptotic behavior during extrapolation. BM2 EOS 
has been experimentally verified for a wide range of minerals under high pres- 


sure. 
(AY om 


e po is a reference density at ambient conditions 


Pgm2(p) = 


e Ko is a reference isothermal bulk modulus at ambient conditions 
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Both po and Ko are inputs in BM2 EOS for a particular crystal phase. If this 
particular phase does not exist at ambient conditions (ordinary temperature 
and pressure), and then appropriate nominal extrapolated values are adopted, 
such as the case for Ice VII. 

Note that the pressure increment dP in BM2 EOS can be written as: 


wesw (2-4 (2) GE) om 


Thus, the bulk modulus K at any density can be written as: 


_ dP 3 7T (pN 5 (DN 
5 mom i ME ; (4) oe 


The bulk modulus increment dK in BM2 EOS can be written as: 


(3) (2) -G) (2) JG) 
dK = - : Ky: . . -d| — 193 
pt © po 3 Po Po eel 
Thus, the derivative of K with respect to P, can be written as: 
E Ge. 
dK 3 Po 3 Po 
K'= = (194) 
06 068" 
3 po 3 po 


The asymptotic behavior of K'(= dK/dP) (dimensionless) is as follows: 


dK 

K aS- =4, (195a) 
dP P=Po0 
dK 

Ke = 7/3 2.33. (195b) 
dP | so 


The behavior of K'(= dK/dP) is illustrated in Fig. 
The initial value of Kj = 4 is close to many experimentally-determined 
values for silicates through static compression. 


4.3.4 Thermal Pressure: 
Debye Model 


On the other hand, we approximate Pihermal With the Debye Model [BG] for T > 
0p, which is generally true for planetary interiors deeper than a few kilometers: 


Eyib 


Pdebye(p; T) =y I (a $ Kr)p : (T —0.23 Op) (196) 


e Op is the Debye temperature 
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Figure 22: K'(= dK/dP) versus p/po 


e a is the thermal expansion coefficient, in unit of inverse temperature 
(KT) 


e Kr is the isothermal bulk modulus, in unit of pressure (GPa) 


where the underseript D of (a-Kr)p means that the product (a:K7) is evaluated 
at the Debye temperature 0p. For T > Øp, (a- Kr) stays approximately 
constant, because of one important thermodynamic identity: 


a:Kr=-p:O0, (197) 
e ~ is the Grüneisen Parameter 
e p is density 
e C, is the isochoric heat capacity 


Above 9p, C, approaches a constant of 3- R per mole of atoms, known as the 
Dulong-Petit limit. This limit is because above the Debye temperature 6p, 
each atom in the crystal lattice realizes its full degree of freedom and becomes 
equivalent to an independent 3-D harmonic oscillator. Therefore, the above 
expression of Pihermal can be further simplified into the following form: 


Pipermail, T) =a: (T = b) (198) 


where a = a: Kr and b = 0.23- Op are constants to be determined, either by 
theoretical approaches or by fit to experimental data. 
In summary, we write the analytic form of pressure in a solid phase as: 


P(p,T) = 5 Ko: H” (Ey ta- (T—b) (199) 
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From this relation of Eq.|199| one can easily calculate the density-pressure (p- 
P) relation for any isotherm (with constant temperature). Because of the nature 
of atomic vibration in a crystal lattice, this equation tells us that regardless of 
the density or pressure, the thermal pressure contribution as described by Debye 
model is only a function of temperature. This allows the separation of variables 
and simplifies our calculations. 


4.4 Ice-X VIII: Superionic Ice 


This sub-section deals with calculating the thermodynamic parameters of supe- 
rionic ice, from the standpoint of BM2 EOS for cold compression, and Debye 
model for thermal pressure, as introduced above. 


4.4.1 Calculating Isentropes of Superionic Ice 


The fit of the Eq. to the experimental and ab initio simulation data of 
superionic ice [64] give the following optimal fit values for the superionic phase: 


Ko = 110 GPa 
po = 2.1 g/cc 
1 GPa 
= 0.0117 GPa/Kelvin = — - — 
a = 0.0117 GPa/Kelvin ZE K 


b = 300 Kelvin 


b 
Op = — 7 1300 Kelvin 


0.23 
K 
Cy = at 4.157 - di (200) 
g 
Since 
dinT dinT 
NE sp PE 201 
AG (sez). p E Peon (201) 
From Eq. 
-K 
Peonst = AT = £ = 2.8145 g/cc (202) 
Gå Cy 


Then, Eq. can be integrated to find the functional relationship between 
T and p along an isentrope: 


T = T) -exp ( = s) (203) 
p 


where Ty is constant of integration and is determined by the initial condition. 
The initial condition refers to the specific entropy s this particular isentrope 
corresponds to. 
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Then, the total pressure along an isentrope is expressed in the following 
analytic form: 


3 7/3 å 5/3 Ge 
Pisentrope(P) = 3 . Ko . ( ) ( ) | | a: (To . e P = b) (204) 


Po Po 


One has to be cautious as this relation of Eq. only applies to within 
a single solid phase. If the isentrope encounters a phase boundary, such as a 
Fluid-Solid phase transition, or a Solid-Solid phase transition, then there will 
be a discontinuity in the specific entropy, and a discontinuity in the density (or 
equivalently specific volume) as well. Then this calculation needs to be carried 
out separately for each phase. 

The isentropes of the superionic ice (Ice-eighteen) can be calculated from 
the ab initio simulation data from [64], as shown in Fig. The local slope of 
isentrope, i.e., the Griineisen parameter y, is calculated from Eq. 

It is also quite interesting to see how the isentropes behave below the density 
of critical point (Fig. (24). This regime people often refer to as Real Gas, as it is 
still considered a gas phase, while the inter-molecular interaction is quite strong 
and cannot be neglected. One such type of EOS that we will introduce later on 
is called the Berthelot EOS. 

Sometimes, it is useful to view these isentropes in a linear-plot of density p 
rather logarithmic. Thus, we come up with the following Fig. 


4.4.2 Calculating Melting Curve of Superionic Ice 


The melting curve here defines the phase boundary between the liquid and 
the superionic ice (also known as Ice-eighteen [30]). It has a shallower slope 
compared to the isentropes in the same regime in the temperature-density (T-p) 
phase diagram. It can be calculated according to Eq. 

With this knowledge in mind, we refine the High-Pressure Phase Diagram 
of H20 with the density jump at the phase transition between the fluid and 
superionic. See Fig. 

These eight isentropes can also be viewed in the p-P (density-pressure) plot 
(Fig. 27), in order to calculated the mass-radius curve corresponding to each 
isentrope. 

The mass-radius curves of these eight different isentropes, as a result, are 
shown in a small-scale mass-radius plot (Fig. and a large-scale mass radius 


plot (Fig. B9). 


4.5 EOS of Ice VI and Ice VII 


According to [25], the EOS parameters for Ice VI and Ice VII are as follows. 
Note that the underscript ”9” suggests extrapolation of the parameter under 
ambient conditions (Pp) =1 bar and Tọ =300 Kelvin). These parameters are 
obtained by fit to BM2 EOS with temperature-corrections as: 


69 


Density-correction: 


WT = po(P, To) - exp ( ape T)) (205) 
Bulk-modulus correction: 
Ok 
K(T) = Ko + (T —T)- (37) (206) 
OT jp 


The melting curve of both Ice VI and Ice VII is fit to a Simon-Glatzel form 


of Equation [73]. 
T C 
Pu = Py +a: (3) -1) (207) 
To 


where a and c are constants to be determined. More recent development 
involves the Kechin Equation [54]. 


4.5.1 Ice VI 
Thermodynamic parameters of Ice VI [25]: 


po 1.27 g/cc 

Ko ~ 14 GPa 

ao © 14.6 - 1075K! 
J/K 

Cy & 2.5 I/R 


-K 
y= VAT 0.65 
Cy * P 


(208) 
Melting Curve of Ice VI [49]: 


T 4.6 
Pu © 0.63 + 0.68 - (( = ) = ) GPa (209) 
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4.5.2 Ice VII 
Thermodynamic parameters of Ice VIT [25]: 


po © 1.45 g/cc 
Ko ~ 20 GPa 
ao ~ 11.6 - 1075K! 
J/K 
Cy & 2 J/K 
g 
a: Kr 


Cy’ P 


= 0.81 


Y == 
(210) 


Melting Curve of Ice VII 1], fit to experimental data up to 60 GPa: 


Tix 4.32 
x 2. : . - —1 P 211 


and the density of Ice VII along the melting curve is [41]: 


Pu Pu 
pm =% 1.45+0.4- ( exp ( 0.0743 ai) ) 428 ( exp ( On) g/cc 


(212) 
Thus, a good rule of thumb for the Griineisen parameter in the regime of 
both Ice VI and Ice VII is that y = 0.7 ~ 0.8. 


4.6 3D EOS of H2O 


Finally, in order to aid our global understanding of the landscape of H20 EOS, 
we visualize the isentropes in a 3D pT P EOS of H20 in Fig. Fig. Fig./32] 
and Fig. 
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Figure 23: Isentropes of fluid-H2O and superionic-H20. The arrows (orange for 
superionic and purple for high-pressure fluid phase) show the local slope of an 
isentrope, that is, the Grüneisen parameter y which is calculated from applying 
Eq. to the ab initio simulation data of [64]. Pressure contours are shown 
and labelled with its value in GPa. The density discontinuity between the fluid 
and superionic phases is not shown in this plot, but will be elaborated in more 
detail in the next section. 
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Figure 24: Isentropes (colored) and Isobars (black) of H2O below and near its 


critical point. 
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Figure 25: Isentropes of H20 in linear-p and In T plot. 
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Figure 26: Eight isentropes of fluid-H2O and superionic-H20, plus the density 
discontinuity at the phase boundary, are chosen here for the purpose of calcu- 
lating mass-radius curves of pure-H2O planets. Pressure contours are shown in 
black. Their smooth behavior suggests a possible approximation by power-law 
for future simplification of calculations. That is y ~ constant. 
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Figure 27: Density-Pressure (p-P) plot of isentropes of H20. The discontinu- 
ity in each curve corresponds to the density jump at the phase boundary be- 
tween fluidic and superionic. The color scheme of the curves are the same 
as Fig. At high-pressure or high-density, all isentropes converge and have 
smooth behavior. This suggests that the EOS in this regime can be approxi- 


mated by I = (gn ) = constant. The achieve higher fidelity, one could add 


s 
a temperature-correction term. 
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Figure 28: Mass-radius curves of eight H20 isentropes below 4 Earth radii (Rg). 
These mass-radius curves suggest that with the appropriate consideration of 
interior temperature/entropy sustained by planet’s own energy source, it could 
help explain some of the exoplanets observed in this regime as water-rich planets. 
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Figure 29: Mass-radius curves of pure H2O isentropes compared with the pure 
Hydrogen isentropes. The color of each mass-radius curve corresponds to differ- 
ent assumed specific entropy. The black mesh provides another set of contours 
of the central densities of planets. 
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Figure 30: A simplified 3-D p — T — P EOS of H20 EOS with isentropes (red- 
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Figure 31: 3-D p — T — P EOS of H20 EOS with isentropes (red-thin curves 
and colored dots) and isobars (black curves) visualized. Recall that along an 


isentrope, I = one , and thus, I are just the tangential slopes of these 3-D 


isentropes when viewed from the perspective of the P-p plane, c.f. Fig. 


Figure 32: 3-D p—T — P EOS of H20 EOS with isentropes (red-thin curves and 
colored dots) and isobars (black curves) visualized. Viewed from the perspective 


of the T-P plane, the tangential slopes of isentropes are % = ee = 


s 
OlnT OlnP 
Olnp Olnp i 
s s 
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Figure 33: 3-D p—T— P EOS of H20 EOS extending down to lower density (p) 
below the critical point (big red dot). This shows how quickly it approaches the 
ideal gas EOS below the critical point: the contours of isochors (red thin curves), 
isobars (black thin curves), and isotherms (green thin curves) all become linear 
in this low density regime in the lg-lg-lg plot. The gas-liquid equilibrium (co- 
existence) area below the critical density (pe) is colored purple. This observation 
reinforces the idea to use a modified ideal gas EOS to describe the entire fluid 
regime of H20. 
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5 EOS of NHs, No, CHyg.,... 
5.1 Ammonia (NH3) 


According to our earlier discussions, the Temperature-Density (T-p) phase di- 
agram of Ammonia should resemble that of water. To the first-order approxi- 
mation, we just need to re-scale the H20 phase diagram properly to obtain the 
one for Ammonia. 

Then, we will implement the ab initio calculation results from Rostock 
Group, in particular, Mandy Bethkenhagen [24] [23]. Their calculations show 
that NH3 also contains a super-ionic phase. 


5.1.1 Critical Point (C.P) for Ammonia 


The analytical formulae of NH3 are anchored at its critical point (c) [9]: 


Pe = 0.226 g/cc 

T, = 405.5 K = 132.3°C 

P, = 1.128- 10” Pa = 112.8 bar 
ne = 0.226/17 = 0.013 mol/ee 


(213) 
Here ne is the molar density of ammonia (NH3) at its critical point. 
5.1.2 Gas-Liquid-Ice Triple Point (t1) for Ammonia 
The Gas-Liquid-Ice triple point (t1) has the following properties [9] [5]: 
pu = 0.733 g/cc for liquid 
Ta = 195.4 K = —77.75°C 
Pri = 6060 Pa ~ 0.06 bar 
ni = 0.733/17 = 0.043 mol/cc for liquid 
(214) 
Here nı is the molar density of ammonia (NH3) at its triple point. 
5.1.3 Re-Scaling 
Now, let’s look at the ratio of densities at the C.P.: 
c NH: 0.226 
PeNHs L lA w 0.702 (215) 


Pemo 0.322 


and the ratio of densities at the Gas-Liquid-Ice Triple Point (t1): 
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Pt1, NH3 — ao x 0.733 


Pt1,H20 1. 


These two ratios are close to one another. 


5.2 Nitrogen (N2) 
5.2.1 Critical Point (C.P) for Nitrogen 


The analytical formulae of Nə are anchored at its critical point (c) : 


pe = 0.314 g/cc 
Te = 126.19 K = —146.96°C 


P. =3.4- 10° Pa = 34 bar 
Ne ~ 0.0112 mol/ee 


Here n, is the molar density of nitrogen (N2) at its critical point. 


5.2.2 Gas-Liquid-Ice Triple Point (t1) for Nitrogen 
The Gas-Liquid-Ice triple point (t1) has the following properties [8]: 


pu = 0.867 g/cc for liquid 
T, = 63.14 K = —210.01°C 
Py = 1.2510" Pa = 0.125 bar 


ni ~ 0.031 mol/ee 


Here ny is the molar density of nitrogen (N2) at its triple point. 


5.2.3 Re-Scaling 
Now, let’s look at the ratio of densities at the C.P.: 


Pena 0.314 
Na NG 
Pemo 0.322 


and the ratio of densities at the Gas-Liquid-Ice Triple Point (t1): 


PN: 4 oar ~ 0.867 


Pt1,H2O 1. 


These two ratios are close to one another. 
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(216) 


(217) 


(218) 


(219) 


(220) 


5.3 Methane (CH,) 
5.3.1 Critical Point (C.P) for Methane 


The analytical formulae of CH, are anchored at its critical point (c) I]: 


pe © 0.0101 - 16 ~ 0.162 g/cc 
T, = 190.6 K = —82.55°C 

P, = 4.61 - 10° Pa = 46.1 bar 
Ne ~ 0.0101 mol/cc 


Here ne is the molar density of methane (CH4) at its critical point. 


5.3.2 Gas-Liquid-Ice Triple Point (t1) for Methane 
The Gas-Liquid-Ice triple point (t1) has the following properties [7]: 


pu 0.439 g/cc for liquid 

Ti = 90.67 K = —182.48*C 

P, = 1.169 - 10? Pa = 0.1169 bar 
ni ~ 0.439/16 ~ 0.0275 mol/cc 


Here ny is the molar density of methane (CH4) at its triple point. 


5.3.3 Re-Scaling 
Now, let’s look at the ratio of densities at the C.P.: 


: 162 
Patia _ 0-162 _ 4503 
Pc,H2O 0.322 


and the ratio of densities at the Gas-Liquid-Ice Triple Point (t1): 


4 
PLO _ 0.489 9.439 
Pt1,H2O 1.0 


These two ratios are close to one another. 


5.4 Discussion 


(221) 


(222) 


(223) 


(224) 


The comparison among all these gas species reveal their similarities in terms of 
their molar densities at their critical points and their triple points. Generally 


speaking, we have for all of them (H2, H20, NH3, No, CH4,-..): 
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Ne ~ 0.01 mol/cc 
ni ~ 0.03 mol/cc 
nu/ne* 3 


(225) 


This approximation is generally good to within a factor of 2. Even more so 
for the ratio of ny /ne- 

This supports the idea of the principle of corresponding states [53], and will 
aid our future discussion of gaseous envelopes. We can think of ne ~0.01 mol/cc 
as a general critical molar density at which cohesion of molecules transition from 
being close to zero to non-zero as pointed out by D. Mendeleev (1870). 

nz, can be understood as another critical molar density at which cohesion 
of molecules has exerted its full power, and any compression beyond this point 
is strongly resisted by the repulsive interactions of molecules, that is, the Pauli 
repulsion at shorter ranges due to overlapping electron orbitals. 

To go from ne to nz, only requires a three-fold compression. That means a 
mere a factor of 3 = 1.44 reduction in the linear dimension or linear spacing 
between neighboring molecules brings about the overturn in the inter-molecular 
interaction from cohesion to (strong) repulsion. 


5.5 An Approximate EOS 


Let’s define the density of liquid at ambient conditions as pọ ~ 1 g/cc. In 
contrast, the density at critical point is typically one-third of that. 

The density under other conditions as p. Their ratio can be defined as a 
dimensionless parameter: 


n = p/po (226) 


Furthermore, we define another dimensionless parameter for temperature: 


T=T/T. (227) 
Then, approximately the pressure can be expressed as: 


2 
ae n 4 ron ep (42) (228) 
This EOS captures both the increase of pressure due to electronic Coulomb 
repulsion, and critical behavior, as well as thermal effects (including both ideal 
gas behavior at low-density high-temperature or atomic vibrations in crystal 
lattices at high-density). 
The reason it works is because at low-density: P = n- kp :T, where n is the 
number density of molecules, where each molecule contains several (~ 3) atoms. 
At high-density: Pinermal © 3-A- kpg-T, where ni is the number of atoms. It can 
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be tuned to higher-accuracy by adding one or two extra parameters inside the 
exponential term. 

Ea. can reproduce the approximate behavior of H20, NH3, N2, CH3 
from ~one-tenth g/cc up to ~5 g/cc, which covers the transition regime of fluid 
from gas-like behavior at lower density to more crystal-like behavior at higher 
density. 


10 
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0.01 | 
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Figure 34: Pressure contours (labelled in GPa), calculated from this approxi- 
mate EOS (Eq. [228). 


87 


5.6 Approximate gas-liquid and gas-solid equilibrium 
and its relation to atmospheric escape 


On the other hand, the pressure of gas-liquid and gas-solid equilibrium can be 
crudely represented as: 


L 
Peq ~ Poo + exp (- 7) (229) 


This equilibrium can be thought of an escape process, similar to the es- 
cape process occurring in the planetary atmosphere. L is the energy that each 
molecule has to overcome in order to go from the condensed phase (solid or 
liquid) into the dilute phase (gas) which is more or less close to ideal gas con- 
ditions. As a crude approximation, we treat L as the same for gas-liquid and 
gas-solid equilibrium. 

If the system is close to thermal equilibrium, then the probability distri- 
bution of molecules (in both the condensed and the gas phases) satisfy the 
Boltzmann distribution. The probability of finding the molecules within the 
energy range of E ~ E + dE would be: 


f(E) -dE x exp (- =) -dE (230) 


The probability of finding the molecules above a certain energy threshold L 
would be: 


a f(E)dE x exp ( = 7) (231) 


In fact, we realize this is exactly the scaling factor multiplying to pressure 
in Eq.|229 

On one hand, with regard to phase equilibrium, L is interpreted as the (nom- 
inal) latent heat of phase transition. P is understood as an extrapolated pres- 
sure as T — oo. For simplicity, L/kp defines a new Temperature Ty = L/R. 
Then, Eq. [229]can be re-written as: 


To 
Pog ~ Po: — — 232 
ä exp ( F) (232) 
For example, for H20, 


Paa ~ 1 TPa 
To ~ 6000 K 


Interestingly enough, the critical point often occurs around To/T ~ 10 and 
the triple point often occurs around To/T = 20. Thus, for H20, 
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T, ~ 600 K 
Ta ~ 300 K 


From the magnitude of Peq ~ 1 TPa, it indicates that it arises from electron 
degeneracy pressure, but as a manifestation of cohesion here. Thus, the factors 
of 10 and 20 within the exponential seem magical! 


Poo | Peritical ~ exp (10) x 2.2 xX 104 
Pes | Pitot ~ exp (20) = 0.5 x 10° 


On the other hand, with regard to atmospheric escape, L can be understood 
as the depth of gravitational potential well at planet surface: 


G: Mpu 


L 
Rp 


(233) 


where p is the mass of each molecule or particle participating in the escaping 
process. We know that escape velocity Vesce is defined as: 


2.G-Mp 
esc = 2 4 
å ( Rp ) ( : ) 


Therefore, L can be re-written as: 
L=1/2-p" Våge (235) 


and, L/kgT is again dimensionless. 
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6 Shock-wave Experiments on Ammonia and Methane: 
What to be expected 


We have put efforts into understanding the Temperature-Density Phase Dia- 
gram of ammonia and methane, in order to guide the shock wave experiments 
of these two species to be carried out on Z-machine at Sandia National Lab- 
oratories. These two chemical species are of particular interest in the current 
age of exoplanetary science, as they together with water may comprise a signif- 
icant fraction of planetary interiors among exoplanets [87]. Therefore, their 
physical and chemical properties in the high-pressure and high-temperature are 
interesting topic to pursue. 

There are existing ab initio calculations of the EOS of ammonia [59] [24], an 
analytic EOS of Ammonia up to ~3000 K and ~20 GPa [50], and existing exper- 
imental data [70], to be compared with, and guide the up-coming experiments. 
See also Fig. 

The Dissociation, Ionization, Super-Ionization, and Metallization of ammo- 
nia are key features to be expected in its phase diagrams. In particular, the 
super-ionic phase has already been experimentally verified for water [64] [63]. 
The super-ionic phase of ammonia or water-ammonia mixture is expected to 
exist [52]. 

In order to hit this super-ionic phase of ammonia in terms of the shock-wave 
experiments, we need to either: 


e use reverberating shocks technique starting with the liquid density at am- 
bient conditions to approximate an isentrope, or, 


e use a single shock starting with a pre-compressed high-pressure ice phase, 
which is already 60% or denser than the liquid 


For methane, a simple analytic EOS was put forth by G. Kerley that goes 
up to ~10* K and ~2.5 g/cc [55], which agrees with the experimental data from 
W. Nellis. A detailed experimentally-verified methane EOS table for relatively 
low-temperature (up to 625 K) and low-pressure (up to 1 GPa) also exists [72]. 

Recent ab initio simulations suggest that methane is expected to dissociate 
first into Hydro-carbons (above ~ 100 GPa), and further into Carbon (Diamond) 
and Hydrogen (above ~ 300 GPa) [73] [42] [£9] M4. 

For future purpose, in terms of their mixtures, ab initio calculations show 
that a linear mixing model for water-ammonia system is good enough approxi- 
mation for a wide range of pressure-temperature phase-space [22], because these 
two species are very soluble in each other, and behave chemically similar even 
under high-pressure conditions. However, ammonia ice and water ice may even- 
tually separate above ~ 500 GPa [71]. 

In terms of the water-ammonia-methane triple mixture, the inclusion of 
methane (or carbon) into the water-ammonia mixture is more tricky at high- 
pressure high-temperature conditions. Methane/Carbon may separate out from 
the mixture, to either rain down as diamonds, or may float upwards as gas as 
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Figure 35: Schematic phase diagram for water and ammonia, or their mixture. 
Also shown are trajectories of isotherm, isentropes, and single/multiple-shock 
Hugoniot. A single-shock Hugoniot will likely miss the super-ionic regime. The 
detailed differences among the various ice structures/phases have been ignored 
as they are all approaching close-packing at sufficient pressures. 


it is somewhat lighter and more volatile. Some ab initio simulations are also 
available [62]. 
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7 Light Gaseous Envelopes 


7.1 Envelope-Core Radius Relation: 
Introducing fugacity f 


This section deals with planets covered with a not-so-massive gaseous envelope, 
and its effect on mass-radius relation. For a light (not-so-massive) envelope, 
typically less than about 20 percent (< 20%) the planet core mass M, (defined 
as whatever underneath the envelope), the calculation is simplified with the 
following assumptions. 


1. First, the gravitational field within the envelope depends only on the radial 
distance r to the planet center, rather than any other conditions (such 
as mass distribution or scale height variation) within the envelope itself, 
because the planet core mass dominates. 


2. Second, the energy transport process within the envelope is assumed to be 
understood and well-defined. There are three modes of energy transport 
process: radiative, convective, and conductive. They together define a 
unique path or trajectory that the thermal dynamical parameters T-p-P- 
entropy... corresponds to one another and to the envelope metric. It is 
exactly along this path that any differentiation is regarded in the following 
differential equations. It is implicitly path-independent. 


Under these two simplifications, the equation governing hydrostatic equilib- 
rium is written as: 


dP G- Me 
Fr i ee E (236) 
By re-arrangement of terms, we have: 
dP G- Me 1 
= r= G-aea(2) (237) 
p r r 


The RHS of Eq. can be integrated straightforwardly from the envelope 
top to the envelope bottom. Let's define this integral as parameter zx: 


bottom 
1 1 1 
=f auaf t) <a. ( = ) (238) 
top r Tbottom Ttop 


Here Re is the radius measured at the core surface or the envelope bottom, 
and Rp is the radius defined at the envelope top, also known as the total planet 
radius. We can re-express the relationship between Rp and Re in terms of x as: 
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Ry = Re: m (239) 
2 2 
(G ` Me/ Re) 


The cluster GMe/Re has its physical meaning. It is the square of escape 
velocity at the core surface. It measures how deep inside a gravitational potential 
well does the core surface resides. It has the unit of specific energy (energy per 
unit mass). 

This Envelope-Core Radius-Relation (Eq. |239) helps understand how the 
mass of a planet core Me may affect the extent of its envelope. Let’s perform 
the following thought experiment: If Re stays constant, then, 


e when Me N; (MTR) TEJ ^ Rp Z 


e when Me /, (Gio) TRJ © Ry \ 


Another corollary of Eq. is that if x becomes too large, then the de- 
nominator turns negative and situation becomes nonphysical. In reality, this 
suggests that the envelope becomes unbound. 

For the LHS of Fa. a we need to consider the integral f dP | Without 
any a priori knowledge of the thermal structure of the envelope, it is assumed to 
be at a constant temperature (isothermal) determined only by the host stellar 
radiation, i.e., the equilibrium temperature Teq. This is similar to the situation 
on Earth surface where the heat flux from interior is orders-of-magnitude smaller 
than solar radiation. The mean temperature of fluids on Earth surface, i.e., 
the hydrosphere and atmosphere, is mostly determined by the average solar 
radiation. 

Then the integral f a2 is a corollary function of EOS P = Plp, Teq], specif- 
ically for an isothermal process at constant temperature Teq, with the proper 
upper and lower limits applied. This integral also has the unit of specific energy, 
that is, energy per unit mass. 

Alternatively, the same integral f dP can be derived by equating the chemical 
potential u at any given height within this light isothermal envelope. In this 
manner, this integral is related to fugacity f, introduced by G. Lewis (1901). 


3 T f aiton bottom dP 
z= R- Teq In | fott = I aF (240) 
H ftop top p 


here the molar mass yz is needed to convert energy per mole to energy per 
unit mass. 

The fugacity f is especially useful when considering an envelope consisting 
of multiple species, and chemical reactions occurring among the constituent 
species. When writing out an equilibrium equation for chemical reactions, fu- 
gacity f is used often instead of partial pressure. f can be viewed as the pressure 
a given real system must produce to have the same action as an ideal system, 
given the proper asymptotic limit (constant of integration): limp_,o £ = 1. 
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f is the corrected (effective) pressure which takes into account the attractive 
and repulsive intermolecular interactions. 

It should be emphasized that any f is only specific to a given temperature 
Teq: For a different Teq, f % needs to be re-calculated. 

As regards to the limits of this integral: 


e The lower limit relates to what pressure-level or density-level a transit 
observation defines the rim of a planet. It is typically on the millibars 
(100 Pa) level or a few tens of millibars level, depending on the chemical 
and thermal conditions in the envelope (please see the detailed calculation 
of the slant-path optical depth which takes into account planet surface 
curvature in [34] and its supplement). 


e The upper limit relates to the total mass or mass fraction of the envelope. 
This relation is elucidated by the next sub-section. 


Once f dP that is, x is computed from an EOS, then we can use it to com- 
pute Rp according to Eq.|239| To aid our calculation, we define a dimensionless 
parameter y as: 


E bottom qp G- Me G- Me 
e a a ee 


Then the Envelope-Core Radius-Relation (Eq. |239) becomes: 


1 
In many cases, it is more convenient to express planet mass and radius in 
Earth units. Therefore,to aid our calculation, we introduce yet another dimen- 
sionless parameter z: 


7 bottom gp G-Mẹ\ _ G : Ma 
ENN 


The denominator equals a constant in unit of specific energy (erg/g): 


i . 1078 2. . 127 
(2 va _ 6.674 - 107%c/g/s? - 5.972 - 10 ee T 


Re 6.371 - 108cm 


Then the Envelope-Core Radius-Relation becomes: 


m 


Or, more simply if Me and Re are already expressed in Earth units: 


(245) 


R, =R, 
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1 


Me 
s- (/ (5) 
Re 
This relation (Eq. |246) is useful when comparing to contours of Mp/Rp 
(planet mass divided by planet radius, in Earth units) in the mass-radius plot. 


(246) 


Rp = Re- 


7.2 Envelope-Core Pressure Relation 


In order to understand the relation between the mass of the envelope Meny = 
(M, — Me) and the pressure Peny it exerts onto the planet core, We now express 
the hydrostatic equilibrium Eq. in variable mass m instead of radius r: 


dP G-m 1 
dm mri  4n-G-m- (gm) eee) 


Notice that here gravitational acceleration g is a function dependent on m or 
r. Its square occurs on the RHS of the above Eq. Through manipulation, 
we have: 


dP =- . (248) 


Eq. [248] tells that Pressure P and Gravity g are essentially two sides of the 
same coin. Pressure arises from gravity from matter itself to compress itself, 
and thus, the square in the equation. 

It is very tempting to integrate Eq. 


g3 M, 
Poy x g í env 
47G Mp 


(249) 


where g? is some kind of mass-averaged value within the envelope. It would 
be nice if we can somehow evaluate g or g? within the envelope. In fact, this is 
feasible, as long as: 


1. the envelope is light (not-so-massive) (< 20% Me), compared to the planet 
core mass Me or the total planet mass Mp = (Me + Menv). 


2. the envelope mass is concentrated towards its bottom, so that gravity g 
can be approximated by its value evaluated at the bottom of the envelope, 
or equivalently, on top of the planet core: ge. 


If both (1) and (2) are satisfied, then the integral from the envelope top 


to the envelope bottom gives an approximate pressure Pny that the envelope 
exerts onto the planet core due to gravitation: 
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and Me + Meny is no different than the total mass of the planet Mp. gc is of 
course related to the core mass Me and core radius Re as: 


_ GM 
%= TR)? 


P E G. M2 Menv =P Menv (252) 
v C \ Aar Ri M. } ° \ Me 


Here we introduce P, as a characteristic pressure in the deep interior of the 
planet core. For example, for Earth, P. = {x 100 GPa or 1 megabar, which 
happens to be the pressure at the core-mantle boundary inside Earth, and about 
one-third of the pressure at the Earth’s center. This is a fundamental relation 
that relates the envelope pressure at its bottom to the characteristic pressure 
deep within the core, regardless of the detailed chemistry and physics of the 
envelope or the core themselves: 


Pi AJ Meny 
()-(4) av 
7.3 Real Gas EOS: 


Introducing Compression Factor Z 


(251) 


Therefore, 


The complexity of the fluid H20 EOS benefits from the multitude of experimen- 
tal data available and its application in industry. The accuracy achieved there 
is typically less than 0.1% level. For other gases, let’s start from scratch and 
solve a simpler problem first. As for our purpose in planetary science, a 10% 
accuracy is tolerable, so let’s construct a EOS from a simpler standpoint. 

To aid our discussion, let’s first define a dimensionless parameter: Compres- 
sion Factor Z as: 


_ P: Vm a P 
© RT (p/u) RT 
Where Vm is the molar volume, and p is the molar weight of molecules. For 
Ideal Gas, Z = 1. Deviations of the Compression Factor Z from unity are due 


to attractive and repulsive intermolecular interactions. By its definition, it is 
linked to the Fugacity f as: 


(254) 
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f= rer ft ar) (255) 


integral performed along an isotherm 


When Z = 1, Eq. |255| immediately gives: 
f=P (256) 


7.3.1 Virial EOS 


Z can also be derived from the first principle through statistical mechanics, 
known as the Virial Equation [53]: 


2 3 
z=1+n (2) +e (2) +p. (2) +a. (257) 
pi m u 


where B, C, and D, ... are the second, third, fourth, etc. virial coefficients. 
These coefficients depend only on temperature. They account for interactions 
between successively larger groups of molecules. They reflect double (B), triple 
(C), quadruple (D), etc. interactions of molecules. The detailed calculations 
are based upon equations establishing the relationships between the energy of 
interaction of molecules, and distance between them, and their arrangement. 

For example, for spherical molecules, coefficient B can be calculated as: 


B= -Na L ( exp ( =) . (47r?) - dr (258) 


where E, (r) is the intermolecular-force potential energy expressed as a func- 
tion of distance r. Na is the Avogadro number and kg is the Boltzmann con- 
stant. B has the unit of volume per mole. One example of E,(r)-potential to 
plug in is the Lenard-Jones 6-12 potential. 

For comparatively low pressures (in terms of the reduced pressure 7 = 
(P/P.) £ 0.5), we can limit ourselves to the second virial coefficient B, which 
takes into account only double collisions: 


= ees 
Z=1+B (2) (259) 


If we include not only B, but also C, then the cubic virial EOS gives sat- 
isfactory results up to reduced density ô = (p/p.) S 2 and inverse reduced 
temperature T = (T/T) < 2: 


2-14 (2) +o. (2) (260) 
u u 
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The cubic virial EOS (Eq. has all the nice characteristics of Van der 
Waals EOS, but without its fatal singularity. The cubic virial EOS accurately 
represents the gas-liquid equilibrium of most substance from the critical point 
down to the triple point, where solid phase starts to appear. 


7.3.2 Berthelot EOS: 
For comparatively low pressures (<100 bar) 
or molar density up to ~0.01 mol/cc 


D. Berthelot (1903) proposed the following EOS [53]: 


Q- T-T 
Z=1 «|1-—6-7? 261 
+ 738 ( å oe 
Here m is the reduced pressure 7 = È, and T is the inverse reduced tem- 


Li 
perature T = 7#. 


7.4 Procedure to Calculate Mass-Radius Relation 


When we calculate the mass-radius curves for a fixed envelope-to-core mass 


Pony 


, we are essentially fixing (3). However, recall the definition 


+, | Menv 
ratio: M. 


of P. that it depends on both Me and Re, and thus, it varies along the mass- 
radius relation (M,-R,) of a given core composition. 

Therefore, we need to follow this procedure to calculate the mass-radius 
relation for a planet core with an envelope: 


e First, we need to identify the mass-radius relation (Me-Re) of the core 
itself. 


e Second, we need to identify the mass Meny or mass fraction Meny /Me of 
the envelope. Since it is a light envelope, we assume that it does not 
exert too much pressure onto the core, and change the core’s mass-radius 
relation (Me-Re). 


e Third, we compute the pressure at the bottom of the light envelope Ppottom: 
G-M? Meny 
ottom = < : 262 
Prott å) (Az) (262) 


e Fourth, carry out the integral f dP from top to bottom of the envelope 
with the appropriate EOS, to find the dimensionless parameter z: 


E Prottom dP Gs Me O G. Me 
(f EJ /( a )--/(G* (263) 
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e Fifth, we apply the Envelope-Core Radius-Relation, to calculate Rp: 


(264) 


R, =R 


Here we show an example of adding a 2-percent-by-mass pure-H»-envelope 
onto a half-ice-half-Earth-like-rocky core and its effect on the mass-radius curves. 
We compare this set of isothermal-envelope mass-radius curves to the Kepler-36 
system and K2-36 system (Fig. (36). 

A minimum point exists as a common feature for these mass-radius curves 
of fixed envelope-to-core mass ratio. It is due to the trade-off between the 
extensiveness of the envelope and the size of the core. Since in the mass-radius 
plot shown here (Fig. [36), mass coordinate is in logarithmic scale, it suggests 
that the overall planet radius R, stays roughly constant over a significant mass- 
range. 

Using this minimum point, we could separate the envelope mass-radius curve 
into a left branch and a right branch. 

Within the left branch, the planet radius R, increases with decreasing core 
mass Me or planet mass Mp. It quickly crosses over contours of constant sur- 
face escape velocity or constant surface gravitational potential (M,/R,=const), 
contours of constant surface gravity (M,/R2=const), and contours of constant 
mean planet density (M,/ R3=const), in the decreasing way. So it will quickly 
become unbound to the planet core and escape. 

Within the right branch, the planet radius R, increases with increasing core 
mass Me or planet mass Mp. In the meantime, the envelope thickness (Rp — Re) 
decreases and then tends towards constant with increasing core mass Me or 
planet mass M,. The planet core holds the envelope tighter at higher mass. 
Thus, the envelope becomes more stably bound to the planet core, and less 
susceptible to escape. 


7.5 Example: Separation of Variables 


Here we give an example to calculate the parameters x, y, and z for a specific 
EOS-a Modified Ideal Gas EOS, and explain the implications of the results. 

For a generic fluid, including both gas and liquid, its pressure-EOS can often 
be treated with the technique of separation of variables: 


P= Prhermal F Presidual (265) 
SVS I 
Part I Part II 
pp Pihermal (T, p) + Presidual (p) (266) 
———— sr 
Part I Part II 


Part I is, for an example, the Ideal Gas EOS. Part II results from inter- 
action among molecules through effective pair potential. High-pressure shock 
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K2-36 b,c (RV): [Fe/H]=-0.15+0:04, age=1.1702Gyr (Damasso 18) 
P „K-36 b,c (TTV): [Fe/H]=-0.20+9:28, age=6.8+1:2Gyr (Carter_1 2) È 
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Figure 36: The K-36 two-planet system and K2-36 two-planet system happen to 
resemble each other by mere coincidence of their designation. The two planetary 
systems, each has one rocky-super-Earth, and another planet with sufficient 
amount of volatiles. Their formation pathways may be similar [33] [27]. The 
planet occurrence rate with Gaia-correction is taken from [27]. 
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wave experiments show that this effective pair potential is quite universal for 
fluid diatomic molecules such as Nə, CO, O» and Hg. It is independent of initial 
fluid densities at ambient pressure and potential parameters that scale with the 
critical points (pe, Te, P) of various fluids [68]. 

Thus, it motivates us to re-express pressure P in terms of molecular molar 
density n (in mol/cc) instead of mass density p (in g/cc), in order to make our 
derivations more generic. This is especially useful when treating a fluid which 
is a mixture of multiple components rather than a single component. In order 
to convert p into n and vice versa, we need to introduce the mean molecular 
weight fi as: 


n = 


(267) 


Els 


Suppose a generic fluid consists of Hydrogen Molecules Hz, Helium Atoms 
He, and other heavier species collectively called Z. Assume it is a Linear Mix- 
ture which obeys volume additive law, then the molar density n follows the 
additive law as follows: 


Ntotal = NH, + NHe + Nz (268) 


Eq. is of course violated if chemical reactions occur within or among 
species which result in association or dissociation (such as Hə dissociates into 
two H-atoms) which change the total number of molecules. However, for now 
we assume this is not the case. This is approximately true for H2 until pressures 
of the order of ~100 GPa, for high pressure will prohibit the dissociation of Hə 
when temperature T is not too high. Furthermore, we assume each molecule 
behaves identically and uniformly in this mixture, and thus, their interaction 
can be described by a single effective pair potential. 

Ln, =2 g/mol, and une= 4 g/mol. uz is the mean molecular weight of the 
heavier species themselves, which is typically on the order of 20 or 30. 


Oe ny, + molar density of Hz (moles per unit volume) 

HH» 

LEE nye + molar density of He (moles per unit volume) 

HHe 
ur nz, < molar density of heavier species (moles per unit volume) 
HZ 


(269) 
However, in astronomy, people often like to component of a gas mixture in 


mass fraction instead of number fraction. The relationship can be expressed as 
follows: 
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X:(mass fraction of Ha) = PH = HH, ` NVH» 


PH2 T PHe + PZ UH» * NH + UHe * Nye + UZ * NZ 
Y:(mass fraction of He) = PHe = HH, ' MH, 
PHe + PHe + PZ UH. ` MH, + UHe * NHe + UZ Nz 
Z:(mass fraction of Z) = pz = HH» ` NH3 
PZ + PHe + PZ MH, ` NH, + HHe + NHe + UZ * NZ 


(270) 
From Eq. the relation between fi and wH,, HHe, Hz can be written as: 


1 X Y Z 
= 4 —+ (271) 
H PH,  HHe Hz 
Now, let’s try the following generic form of pressure-EOS as motivated by 
real experimental data: 


ï 
Pa fi | ——— (272) 
— 1 mol/cc 
Part I 
Part II 


Å is a constant in unit of pressure. It is about 10? GPa or 1 TPa, or 
10!3dyn/cm? in cgs units, for H> and other molecules alike approximately. Ac- 
cording to our demonstration earlier, we treat I' as constant, which is ~2 for 
metallic fluid H2, and H2-He Cosmic Mixture (X = 3, Y= 1), and other 
molecules alike approximately [79]. At low number density n, the temperature- 
dependent Part I dominates. At high number density n, the temperature- 
independent Part II dominates. Then, the integral f a for an isothermal 
envelope at constant temperature Teq and constant mean molecular weight fi 
can be performed: 


bottom bottom 
dP 1 dP 
de | eee | ue (273) 
top P H top n 
Then, in cgs units, 
dP=R-T:dn+à- T-n™t!.-dn (274) 
Sapa NTE 
Part I Part II 
Then, 
dP d 
OG Pn! can (275) 
nn) 
ý EE Part II 
Part I 
Then, 
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bottom bottom bottom 
dP d 
Í Zorr f Saar f n”? . dn (276) 
t t n t 


op É op op 


Part I Part II 
Then, 


bottom 
dP ottom r — 
i — =R-T-In (2e ) HF EN = nn (277) 
ee 


op 
— esl‘ 
Part I Part II 


Then, Parameter x, which is a represent of envelope thickness, is: 


R Nb Å Tr 
q ottom T-1 r-1 
g= — - ln | —— += SA =n (278) 
= = bottom to 
H Ntop H P 7 
————— 
thermal contribution to envelope thickness intrinsic thickness of envelope 


Or, equivalently, since Nbottom > Ntop, 


h- Te ottom A r SS 
mafaa) aapi (Ei) a 


thermal contribution to envelope thickness intrinsic thickness of envelope 


a 
Q 


The interpretation of this result is that Part I gives rise to the thermal 
thickness of an envelope, which is characterized by temperature-dependent scale- 
height. Part II gives rise to the intrinsic thickness of an envelope, regardless of 
the thermal conditions and energy transport process. This result is very general. 

This Separation-of-Variable treatment allows us to understand how a fluidic 
envelope expands under increased stellar insolations, applicable to many short- 
period exoplanets, such as hot-Jupiters and hot-Saturns. 


R e Te ottom A r = 
LE EE (280) 
H Ntop H Ped 
eae S mm 
Part I Part II 


7.5.1 Discussion 


Now, let’s plug in the numbers and compare the actual results, in cgs units 


(erg/g): 
. 101 i Teq vi Nbottom +2. 1013 x Nbottom Bi / (281) 
102K Ntop mol/cc 8/8 
$< 


mm aal N 
Part I Part II 


Ble 
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Ntop corresponds to the number density of the top of the envelope that 
the transit observation probes. Let's treat it as 1073 mol/cc. The exact 
value shall depend on the photo-chemistry of the envelope top. This typically 
corresponds to pressures of millibar level (please see the detailed calculation of 
the slant-path optical depth which takes into account planet surface curvature 
in [34] and its supplement). 

Note that Part IT becomes comparable to Part I in magnitude when npottom 
approaches ~ 107? mol/cc or ~0.01 mol/ce. 

Recall that D. Mendeleev (1870) pointed out that it is useful to compare the 
molar densities among different chemical species, not at their boiling points, but 
at their critical points when the cohesion of the molecules transition from being 
close to zero to non-zero there [53]. Indeed, in previous sections we see that: 


Nc Ho & 0.015 mol/ce 
NoHo & 0.018 mol/ee 
Nc, NH; © 0.013 mol/cc 
Ne, Na © 0.0112 mol/ce 
Ne, CH, œ 0.010 mol/cc 
(282) 
Thus, near or above this critical molar density of ~ 0.01 mol/cc, the cohesive 
interactions among molecules start to take over and turn the substance into a 
condensed state. They should be considered liquid or super-critical fluid from 
there on, and they should not be considered gas-like above this critical molar 
density. 
Using Eq. |244| we can obtain the dimensionless parameter y from z: 


1 Te Nbottom Nbottom Me 
~ —-( 0.133 | = | -1 16- L | 4 MaRe 
á Å Gå) a( Trop )- i (zs) | / (36). ft 
es 


Part I Part II 


(283) 
and recall that: 


1 
Rp = Re: Iy (284) 

Now, let’s focus on Part I which varies from ntop ~ 1078 mol/cc to Nbottom ~ 
107? mol/cc. This is approximately one-million-fold contrast in density. 

To visualize this contrast, let’s imagine a small volume box containing a 
certain number of molecules at the very top of the atmosphere defined by the 
transit observation. Now, when we move this box containing the same number 
of molecules down to the bottom, each side (linear dimension) of it shrinks by 
one-hundred fold, and this gives rise to one-million-fold increase in density. 

On the other hand, to go from the bottom of the atmosphere (~ 0.01 mol/cc) 
all the way to the center of a super-Earth exoplanet (~ 0.1 mol/cc) the density 
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ine 


increase is only about one-order-of-magnitude (a factor of 10 or so). This shows 
how different and difficult it is to compress a condensed phase compared to a 
gas phase. 


Since In | 42 | = In (108) = 13.8, 


10-8 
2 Ti Me 
n= (3) (r) / (æ) a= 


From Fig. [36] one can see that most interesting super-Earth exoplanets have 
M,/Rp Z 3. That suggests Me/Re Z 3 for most of these planets. Moreover, 
G > 2 and it takes the minimum value only when it is a pure Ho atmosphere. 
Moreover, most of these interesting super-Earth exoplanets are < 1000 K equi- 
librium temperature. Therefore, 


Te 
yı <0.3- (ek) (286) 


From this analysis, we show that though gas has low density, it is also very 
compressible under a gravitational potential gradient. It does not take a long 
distance for a gas-like species to build up its pressure and molar density until 
reaching its critical molar density, where the molecular cohesion takes over and 
turns it into a condensed phase. 

As a result of that, fundamentally speaking, a gaseous atmosphere can only 
account for a surface-layer of a planetary body, a mere skin-effect. This state- 
ment even applies to a gas giant like Jupiter or Saturn, because most of their 
interior is considered high-density super-critical molecular fluid and then metal- 
lic fluid rather than a gaseous state, in terms of condensed matter science. 

The fractional boost gained from such a gaseous atmosphere is typically on 
the order of 10% to 30% (10-30 percent) of the planet core radius (Re) if the 
atmosphere composition is dominated by H2-He, and is typically on the order 
of 1% to 3% (1-3 percent) of the planet core radius (Re) if the atmosphere com- 
position is dominated by Oxygen-Nitrogen-Carbon-species since their molecular 
weight is about one order of magnitude higher. 

Then in Fig. the group of super-Earth exoplanets with M,/R, 2 3 
and with lower-than-rocky density raise an interesting question. They may be 
explained by: 


e intrinsic lower-than-rocky density in the fluid-layer of these planets, or, 


e intrinsic heat source or heat reservoir that maintains a high-entropy state 
in the interior of these planets. 


Last but not least, due to proximity to host star, the stability of a such 
gaseous atmosphere on a short-period exoplanet needs to be examined carefully. 
Very plausibly, if such a gaseous atmosphere exists, then it is experiencing escape 
and must be continuously replenished by a reservoir such as the fluid-layer 
underneath. 
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8 Conclusion 


I always thought something was fundamentally uniform with the universe. [82] 
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